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Abstract 

We present the Monte Carlo event generator /C/C version 4.13 for precision pre- 
dictions of the Electroweak Standard Model for the process e + e~ — > ff + nj, 
f = fj,,r, d, u, s, c, b, at centre-of-mass energies from r lepton threshold to 1 TeV, 
that is for LEP, SLC, future Linear Colliders, b, c, r-factories, etc. Effects due to 
photon emission from initial beams and outgoing fermions are calculated in QED 
up to second order, including all interference effects, within Coherent Exclusive Ex- 
ponentiation (CEEX), which is based on Yennie-Frautschi-Suura exponentiation. 
Electroweak corrections are included in first order, with higher-order extensions, 
using the DIZET 6.21 library. Final-state quarks hadronize according to the parton 
shower model using JETSET. Beams can be polarized longitudinally and trans- 
versely. Decay of the r leptons is simulated using the TAUOLA library, taking into 
account spin polarization effects as well. In particular the complete spin correlations 
density matrix of the initial-state beams and final state r's is incorporated in an 
exact manner. Effects due to beamstrahlung are simulated in a realistic way. The 
main improvements with respect to KORALZ are: (a) inclusion of the initial-final 
state QED interference, (b) inclusion of the exact matrix element for two photons, 
and (c) inclusion of the transverse spin correlations in r decays (as in KORALB). 
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PROGRAM SUMMARY 



Title of the program: /C/C version 4.13. 
Computer: 

any computer with the FORTRAN 77 compiler and the UNIX operating system 
Operating system: 

UNIX, program was tested under AIX 4.x, HP-UX 10.x and Linux 

Programming language used: FORTRAN 77 with popular extensions such as long names, 

etc. 

High-speed storage required: < 10 MB 

No. of cards in combined program and test deck: about 21,800, without JETSET, 
TAUOLA and PHOTOS. 

Keywords: Quantum electrodynamics (QED), Standard Model, electroweak interactions, 
heavy boson Z, spin polarization, spin correlations, radiative corrections, initial-state ra- 
diation (ISR), final-state radiation (FSR), QED interference, Monte Carlo (MC) simula- 
tion and generation, coherent exclusive exponentiation (CEEX), Yennie-Frautschi-Suura 
(YFS) exponentiation, LEP2, linear collider, TESLA. 

Nature of the physical problem: The fermion pair production is and will be used as an 
important data point for precise tests of the standard electroweak theory at LEP and 
future linear colliders at higher energies. QED corrections to fermion pair production 
(especially r leptons) at c-quark and 6-quark factories has to be known to second order, 
including spin polarization effects. The Standard Model predictions at the per mille 
precision level, taking into account multiple emission of photons for realistic experimental 
acceptance, can only be obtained using a Monte Carlo event generator. 
Method of solution: The Monte Carlo methods are used to simulate most of the two- 
fermion final-state processes in e + e~ collisions in the presence of multiphoton initial-state 
radiation. The latter is described in the framework of exclusive coherent exponentiation 
(CEEX) based on Yennie-Frautschi-Suura exclusive exponentiation (YFS/EEX). CEEX 
treats correctly to infinite order not only infrared cancellations but also QED interfer- 
ences and narrow resonances. The matrix element according to standard YFS exponen- 
tiation is also provided for tests. For quarks and r leptons, the appropriate simulation of 
hadronization or decay is included. Beam polarization and spin effects, both longitudinal 
and transverse, in tau decays are properly taken into account. 

Restrictions on the complexity of the problem: In the present version, electron (Bhabha), 
neutrino and top quark final states are not included (they will be in a future version). 
Additional fermion pair production is not included. Third-order QED corrections in 
leading- logarithmic approximation are included only in the auxiliary YFS/EEX matrix 
element (which can be activated with the help of input parameters). Electroweak correc- 
tions should not be trusted above the t-quark threshold. The total cross section for light 
quarks for y/s < 10 GeV requires an improvement using experimental data. 
Typical running time: On the IBM PowerPC M43P240 installation (266 MHz, 65 CERN 
units) 4 sec per constant- weight event are needed. This result is for a default/recommended 
setting of the input parameters, with all hadronization/decay libraries switched ON. 
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1 Introduction 



Monte Carlo (MC) event generators have the double purpose to compensate for detec- 
tor inefficiencies and to provide theoretical predictions for distributions and integrated 
cross sections. The second task is more important and difficult. Precision predictions 
of the Standard Model, with a total error below 0.5%, for the process of production of 
the fermion pair in electron-positron scattering (e~e + — > fJ>~fi + ) were first obtained for 
energies close to the Z resonance with the MC event generator KORALZ |IJ (see also 
the latest version ||). The prototype of the modern MC event generator for this process 
was constructed earlier ||, but it could not deliver sub-per cent precision, because it did 
not include electroweak corrections and second order QED corrections (it later became 
a part of the KORALZ package). KORALZ was originally developed for a simulation of 
the r-pair production and decay, and later on was extended to muon, quark and neutrino 
pairs, e~e + — > ff,withf = fi,r,d,u, s,c,b,v. The Bhabha scattering e~e + — > e~e + was 
never included in KORALZ and a dedicated precision MC event generator BHLUMI j|,|5| 
was developed for this process. BHLUMI, at the expense of specializing to small scat- 
tering angles, could deliver at LEP1 energies the integrated cross section with the record 
total precision of 0.06% ||. It should be stressed, however, that the high-precision level 
of KORALZ and BHLUMI was achieved thanks to exclusive exponentiation |7|] (EEX) 
based on the classical work of Yennie-Frautschi-Suura (YFS) ||, in which the multiple 
soft and hard real photons are treated in a completely realistic way, i.e. four-momenta 
are generated, and the infrared (IR) cancellations between real and virtual soft photons 
occur exactly to infinite order. 

At the end of LEP2 operation the total cross section for the process e~e + — > ff will 
have to be calculated with a precision of 0.2%-l%, depending on the event selection |§. 
The arbitrary differential distributions also have to be calculated with the corresponding 
precision. In future linear colliders (LCs) the precision requirement can be substantially 
stronger, especially for the high luminosity option, as in the TESLA case. The above new 
requirements necessitate the development of a new calculational framework for the QED 
corrections and the construction of new dedicated MC programs. The present work is an 
important step in this direction. 

The main limiting factor that prevents us from getting more precise theoretical pre- 
dictions for the e~e + — > ff process is higher-order QED radiative corrections (the QED 
part of electroweak Standard Model). In order to achieve the 0.2% precision tag, the 
virtual corrections have to be calculated up to 2-3 loops and the multiple bremsstrahlung 
up to 2-3 hard photons, integrating exactly the multiphoton phase space for the arbitrary 
event selection (phase-space limits). 

The MC event generators KORALZ and BHLUMI, although representing the state of 
the art of MC evaluation of QED radiative corrections for e~e + — > // at the beginning 
of the LEP2 run, are strongly limited in their development towards higher precision. The 
main limitation is rooted in the use of the spin-summed differential cross sections in the 
otherwise so successful YFS/EEX - instead of using spin amplitudes directly^. For this 

1 This was a sensible choice; for instance, it has saved precious CPU time, which was a big problem 
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reason, certain interferences such as the one between initial-state radiation (ISR) and 
final-state radiation (FSR) had to be neglected, especially in the calculations beyond 
first order. Also, in any processes with more Feynman diagrams, such as the Bhabha 
process, EEX suffers from the proliferation of the interference terms, especially beyond 
first order. EEX neglects ISR-FSR interferences. However, they were often unimportant 
- at the Z peak much below 1% of the integrated cross section. The analogous QED 
interferences among photons emitted from electron and positron in small-angle Bhabha 
were also small, of order 0.01%, at the small-angle range of the luminosity measurement. 
Further improvement on the precision for off-Z-peak (LEP2) and for large-angle Bhabha 
scattering definitely requires reintroduction of these interferences. This is achieved in a 
natural way by reformulating the exponentiation entirely in terms of spin amplitudes. This 
turns out to be quite a non-trivial task, and it was done only recently, see refs. fllOlllll . The 
resulting new "reincarnation" of the YFS exponentiation, called the coherent exclusive 
exponentiation (CEEX), was born. The CEEX scheme is implemented in the present 
/C/C MC program for the first time. In fact we have described in ref. [ID] only the first 



0(a 1 )cEEX version, with the pure QED matrix element, while in the present program we 
have already implemented the bulk of the QED 0(« 2 )ceex matrix element, and also 0(a) 
electroweak (EW) corrections. This new, C(a 2 )cEEX, important development of CEEX, 
together with the wealth of numerical results, will be published separately fl2| . 



in the MC calculations, a decade ago. 



5 



1.1 Ultimate MC event generator for two-fermion final states 

Having briefly introduced the reader to the history and the main characteristics of the 
subject of the MC event generators for the e~e + — > // process, let us come to the 
important question: What is the most complete list of requirements that the precision 
MC event generator for the e~e + — > // process should fulfil in order to satisfy the needs 
of the experiments in the present and future e~e + colliders, that is for the entire centre- 
of-mass energy range from the r production threshold up to 1 TeV? This would cover 
experiments at LEP, SLC, 6-factories, c-factories, r-factories, future linear colliders such 
as TESLA, JLC and NLC, and also high-luminosity experiments at the Z resonance, the 
so-called Z-factories. 

In the following we try to answer the above question, and complete the list of the 
desired features of the ultimate MC event generator for two-fermion final states: 

• The total precision of the integrated cross section has to be at least 0.2%. From 
the above, it automatically follows that we need for the QED part and electroweak 
corrections the entire first-order 0(a) and the QED second-order 0(a 2 ), at least in 
the leading- log (LL) approximation with the exponentiation, e.g. with the 0(a 2 L 2 ) 
contributions, especially for ISR. In fact the QED second-order subleading 0(a 2 L 1 ) 
and third-order leading 0(a 3 L 3 ) are also mandatory, at least for the discussion of 
the theoretical error, but preferably present in the actual MC matrix element. The 
inclusion of fff'f' four-fermion final states will often be necessary, especially the 
production of the additional soft fermion pair of light fermions. The Z-factory option 
with 10 9 statistics would be the most demanding experiment, asking for precision 
better than 0.02%! Note that here we do not attempt to define individual precision 
requirements for all kinds of distributions and averages, such as charge and spin 
asymmetries. We assume that their precision should correspond to a precision of 
0.2%-0.02% for the integrated cross sections. 

• The other role of MC event generators is to provide a realistic picture of the process 
for the detector studies. It means that it is highly undesirable to integrate out certain 
final-state topologies. For instance we know that contributions from collinear real 
photons and from virtual photons combine in such a way that the net effect in the 
integrated cross section is zero or negligible. Nevertheless, for experiments it is of 
vital importance to have all real photons manifestly in the MC event, even very soft 
and very collinear ones: it helps enormously to understand correctly the operation of 
a detector, and hence to reduce the experimental systematic errors. The same is true 
for the emission of additional light-fermion pairs. The above maximal exclusivity 
requirement leads directly to exponentiation in the MC event generator, that is to 
a procedure in which the perturbative expansion is reorganized in such a way that 
contributions from IR real and virtual singularities cancel to infinite order. The 
remaining non-IR corrections are calculated order by order - for example in the 
present /C/C MC the non-IR corrections are included in 0(a 2 ). 

• For the initial beams and outgoing unstable r leptons and t-quarks spin polariza- 
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tions, including all spin correlations, have to be fully taken into account, also in the 
presence of the real bremsstrahlung photons. For beams the longitudinal polariza- 
tion of both beams is the minimum requirement. For decaying final fermions, longi- 
tudinal polarizations are often not sufficient - transverse polarizations are necessary, 
both for r leptons and for t-quarks. In fact not only polarizations are required, but 
also the effects of exact spin correlation between spin polarization vectors of two out- 
going fermions. This again makes it mandatory to the use spin amplitudes and/or 
spin density matrices for the fermion production and decay processes. 

• Non-QED electroweak corrections to the e~e + — > // process are usually calculated 
separately and combined later on with QED corrections in an manner that is, beyond 
0(a), an "ad hoc" recipe. It is not done at the level of the spin amplitudes, but 
rather in terms of the inclusive distributions, having integrated some photon phase 
space e.g. transverse momenta beforehand. Typically such an ad hoc procedure is 
based on the second- or third-order "structure functions" for the incoming electron, 
adding subleading first- and second-order corrections. It can be questioned whether 
such an approach is really justified at the precision level of 0.2%, because it is 
too far from the solid environment of the Lagrangian, Feynman diagrams and the 
exact phase space. The Monte Carlo based on the spin amplitudes offers a natural 
realization of such an environment. It is therefore the only viable solution for the 
problem of combining QED and non-QED corrections beyond first order, without 
any unnecessary ambiguities. Such a Monte Carlo can then be used to cross-check 
(calibrate) the popular semi-analytical approaches, which employ all kinds of ad hoc 
recipes. This important role should not be underestimated, as the semi-analytical 
programs have many advantages of their own and will always be around. 

• Effects due to beamstrahlung will be present at future linear colliders and thus should 
be implemented in the MC event generator. The beamstrahlung structure functions 
should be a "user function", supplied (or replaced) easily by the user, without 
any loss of the efficiency of the MC program. Note that the "luminosity energy 
spectrum" will be known from machine simulation only to a certain extent; the 
true distribution will have to be determined from the inspection of the e~e + — > // 
process, most probably for / = e, at small (~ 1°) or intermediate (~ 10°) angles. In 
such a case, reliable MC predictions of the Standard Model (SM) integrated cross 
sections and distributions in the presence of beamstrahlung, with the non-trivial 
event selection criteria, are of vital importance. In this context, the Z radiative 
return, i.e. the e~e + — > subprocess, may also play an important role. 

• The MC event generator should be maximally upgradable to other processes. Al- 
though, in this specification, we concentrate on the precision SM prediction for the 
two-fermion final states, the MC program should be constructed in such a way that 
inclusion of the other similar standard processes such as e~e + — > W + W~ or non- 
standard processes (supersymmetric) should be relatively easy. Also the change of 
electron beam to muon beam should be possible. This requires in practice that the 
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MC program be constructed right from the beginning in a highly modular way. The 
Object-Oriented Programming (OOP) approach would help. 



• For quark final states, the photon emission should be well combined with the quark- 
gluon parton shower, properly taking QED and QCD corrections at the NLL level. 

• Last but not least the program should also have the option to run in a maximally 
inclusive mode in which the type of the final state (fermion) is chosen randomly, 
event per event, exactly according to its integrated cross section. This requirement 
is not as trivial as it may seem, since the integrated cross section, in the MC, is 
typically known at the end of the run (from the average MC weights). 

Let us note that the above specification goes far beyond the very ambitious (at the time) 
specification of the "ultimate MC" formulated at the end of the 1989 LEP1 Workshop. 
How far are we with /C/C MC on the road to this ultimate goal? 



Feature 


KORALB 


KORALZ 


JCJC now 


JCJC 2000 


QED type 


0(a) 


EEX 


CEEX, EEX 


CEEX, EEX 


CEEX(ISR+FSR) 


none 


none 


{a, ctL; a 2 L 2 , a 2 L 1 } 


{...a 2 L 1 ;a 3 L 3 } 


EEX(ISR*FSR) 


none 


{a, otL, a 2 L 2 } 


{ot, otL, ol 2 L 2 , a 3 L 3 } 


{...a 2 L 2 ,a 3 L 3 } 


ISR-FSR int. 


0(a) 


o(«) 


{ot, aL} CEEX 


{ot, «L} CEEX 


Exact bremss. 




1, 2coll. 7 


1, 2, 3coll. 7 


up to 3 7 


Electroweak 


No Z-res. 


DIZET 6.x 


DIZET 6.x 


New version? 


Beam polar. 


long+trans. 


longit. 


long+trans. 


long+trans. 


r polar. 


long+trans. 


longit. 


long+trans. 


long+trans. 


Hadronization 




JETSET 


JETSET 


PYTHIA 


r decay 


TAUOLA 


TAUOLA 


TAUOLA 


TAUOLA 


Inclusive mode 




No 


Yes 


Yes 


Beamstrahlung 




No 


Yes 


Yes 


Beam spread 




No 


Yes 


Yes 


vv channel 




Yes 


No 


Yes 


ee channel 




No 


No 


Yes 


it channel 




No 


No 


yes? 


WW channel 




No 


No 


yes? 



Table 1: List of features of the present JCJC MC, compared to older MCs for fermion pair 
production, and future plans. 
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1.2 How far are we on the road to the Ultimate MC? 



The present specification of the /C/C MC is summarized in Table [II where we have also 
given those of KORALZ and KORALB for comparison. 

As we see, the present /C/C MC has already all the functionality of KORALB. The 
/C/C MC fulfils completely our ultimate specification for the spin treatment. Up to now, 
the first and the only MC event generator fulfilling the above specification was KORALB 



F5, nfl. However, KORALB is limited to ^fs < 30 GeV, because of the lack of the 
Z-resonance and it does not include more than one photon emission - it is based on 
the pure first-order QED calculation, without exponentiation. In KORALZ, longitudinal 
polarization effects are implemented at the level of the differential distribution for initial 
beams and for outgoing fermions (including longitudinal spin correlations). Effects due to 
transverse spin polarizations are omitted. In the present /C/C program these limitations 
are removed - the complete longitudinal and transverse polarizations are implemented 
for beams and outgoing fermions, including all spin correlations exactly, as in KORALB, 
also in the presence of multiple real photons. 

As compared to KORALB, there is still one improvement to be done in the /C/C 
MC: the CEEX matrix element does not generally rely on the assumption m T <C y/s, 
but certain parts of the actually implemented virtual corrections may still rely on this 
approximation, notably the spin amplitudes for the 7-7 box. This approximation is not 
really necessary and will be corrected in the future version. 

As compared to KORALZ, the /C/C MC covers all its functionality, except for the 
presence of the neutrino channels. We plan to implement the neutrino channel in the next 
version. The most important new features in the present /C/C with respect to KORALZ is 
the inclusion of ISR-FSR interference, which at LEP2 modifies the total cross section and 
charge asymmetry by about 2%, the inclusion of the second-order subleading corrections 
and the inclusion of the exact matrix element for two hard photons. 

The /C/C MC is the first full-scale event generator to include the bulk of second order 
NLL corrections, and it may easily incorporate second-order NNLL corrections. In fact the 
exact double bremsstrahlung and exact two-loop virtual corrections are already included, 
but for the moment there is no need to complete the missing NNLL corrections since they 
are ~ 10" 5 . 
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2 Physics content 



The present version of /C/C MC includes first-order 0(a) QED and electroweak corrections 
and almost complete 0(a 2 ) QED corrections due to the emission of photons from the 
initial- and final-state fermions. It does not include the emission of an additional fermion 
pair. 



2.1 Two types of QED matrix element 

The QED part of the calculation is based on the new 0(a 2 ) calculation with coherent 
exclusive exponentiation (CEEX) at the amplitude level. The older QED matrix element 
based on exclusive exponentiation (EEX) at the differential distribution level (amplitudes 
squared and spin summed) of the BHLUMI type is still present and is used as a backup 
solution, for various tests. In particular the EEX matrix element includes 0(a 3 L 3 ) cor- 
rections, which are still absent in our CEEX; it therefore provides a useful estimate of 
these corrections. 



2.1.1 CEEX spin amplitudes 

Defining the Lorentz-invariant phase space as 



/r n n J 3 

dLlPS n (P;p 1 ,p 2 ,...,p n ) = y(2vr)V 4 )(P-g^n(^^o' W 

we write the 0(a r ) CEEX total cross section for the process 

e~(p a ) + e + ( Pb ) -> f(p c ) + f(p d ) + j{h) + 7 (fc 2 ) + ... + 7%), n = 0, 1, 2, n, (2) 

with polarized beams and decays of unstable final fermions sensitive to fermion spin 
polarizations, following refs. [fL0~l , |12|l , as follows: 

oo 

r (r) = 

flux( 



n l ' ' ~ TX ^ / ^LIPS n+2 (p a + Pb',Pc,Pd, fa,..., k n ) PcEExiP 

a: Pbi Pa Pdi kxi ■ ■ ■ j kn 



(3) 

where 

p { cLx(Pa,P b ,Pc,P d , h, fa, ... , kn) = le^>->^ 0(Q) E 

(4) 



£ 1 £ 2 A a \ a \ b \ b n U<Ti<ra " " " vJ l JJl n Kx<ncr 2 ■ ■ ■ aj J A c A c A d A d 3 

i,j,l,m=0 



in 
4 ) 



where, in order to shorten this and other formulas, we use a compact collective notations 

(P\ _ (PaPbPcPd\ 

\x) — \\ a \ b \ c \ d ) 
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for fermion four-momenta pa and helicities Aa, A = a,b, c, d, 

For fc = 1,2,3, a k are Pauli matrices and a° XfjL = is the unit matrix. The com- 
ponents i{,E2,j,k = 1,2,3, are the components of the conventional spin polarization 
vectors of e~ and e + respectively, defined in the so-called GPS fermion rest frames (see 
ref. Jll] for the exact definition of these frames). We define e Q A = 1 in a non-standard 
way (i.e. pa • £a = me,A = a,b). The polarimeter vectors he are similarly defined in 
the appropriate GPS rest frames of the final unstable fermions (pc ■ he — rrif,C — c,d). 
Note that, in general, he may depend in a non-trivial way on the momenta of all decay 
products, see refs. HTT1, p~5| |17|| for details. We did not introduce polarimeter vectors for 
bremsstrahlung photons, i.e. we take advantage of the fact that all high-energy experi- 
ments are completely blind to photon spin polarizations. See next subsection for more on 
spin effects. 

WE define the complete set of spin amplitudes for emission of n photons in C(a r )cEEX, 
r = 0, 1, 2 as follows: 



3*8° es ■••£)= e n 4r } ^ 0) & > ( 5 ) 

i=i 



n oiX) (P k j. x 



i=l I j=l %] 



m>(2) /pfci fc„\ 

" /Jv n \Xa 1 ■ ■ ■ a n ) 




2-^< Apj} ^ Jpj}Jpi} 



S [j] l<j<l<n S [j] 5 [l] 



(7) 



The coherent sum is taken over set {p} of all 2 n partitions - the partition p is defined 
as a vector (pi, p 2 , . . . , p n ); pi = 1 for an ISR and pj = for an FSR photon; see the 
analogous construction in refs. |I^,|19|. The set of all partitions is explicitly the following: 



M = {(0,0,0,...,0), (1,0,0,. ..,0), (0,1,0,. ..,0), (1,1,0,. ..,0), ... (1,1,1,..., 1)}. 

The s-channel four-momentum in the (possibly) resonant s-channel propagator is X p = 
Pa+Pb~J2i=l Pi K 

At 0(a r ) we have to provide functions flj[\k = 0, l,...,r, from Feynman diagrams, 
which are infrared-finite by construction || . Their actual precise definitions can be found 
in refs. [[TI],[T!|]. Here we shall define only the most essential ingredients. The lowest-order 
/?o°^ are just Born spin amplitudes times a certain kinematical factor 

^ = J^r (8) 
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The Born spin amplitudes 53 ( p x ; X) and other spin amplitudes are calculated using the 



spinor technique of Kleiss and Stirling (KS) [^HH reformulated in ref. |ll[] (GPS). Soft 



factors , u = 0, 1, are complex numbers, see ref. ]TTJ for exact definitions; here we only 
need to know their absolute values 



e 2 Ql 



Pa 



Pb 



e 2 Q 2 



f Vc 



k%p c kiPd 



(9) 



The factor 0(fl) defines the infrared (IR) integration limits for real photons. More 
precisely for a single photon, complementary domain Q includes the IR divergence point 
k = 0, which is excluded from the MC phase space, we define a characteristic function 
Q(Q,k) = 1 for k <E Q and Q(Q,k) = for k ^ Q. The characteristic function for the 
phase space included in the integration is Q(Q,k) = 1 — Q(Q,k). The characteristic 
function for all photons in the M C phase space is 



(10) 



i=i 



In the present program we opt for an Q traditionally defined by the photon energy cut 
condition k° < E m[n . Consequently, the YFS form factor || reads 



Y(tt;p a , ...,p d ) = Q e Y n (p a ,p b ) + Q f Y n (p c ,p d ) 

+ QeQfY n (p a ,Pc) + QeQfY n (pb,Pd) - QeQfY n (p a ,p d ) - QeQfYn(p b ,p c 

where 



Yn(pi,p 2 ) 



2aB{Q, Pl ,p 2 ) +2aUB( Pl , P2 ) 
1 f d 3 k 



= -2a 



P2 

kp2 



(12) 



+ 2aK 



d 4 k 



2pi + k 2p 2 -k 



k 2 (2tt) 3 \2kpi + k 2 2kp 2 -k 2 



is given analytically in terms of logarithms and Spence functions [0. As we see, the 
above YFS form factor includes terms due to the initial-final state interference (IFI). 
In the MC with exponentiation it would be possible to do without Q (declare it empty) 
and rely uniquely on the IR regularization with small photon mass m 7 only |[L0|| . In that 
case formulas ( |T2j ) for the YFS form factor would include only the second virtual photon 
integral. 



2.1.2 EEX differential distributions 

In the case of exclusive exponentiation (EEX), we deal with a spin-summed differential 
distribution. We neglect the IFI and each photon is therefore attributed either to the 
initial or to the final state fermion pair. The initial state photon momenta are denoted 
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by hi, i = 1, 2, n, and final state momenta as i — 1, 2, n' . We shall consider the 
process e~(pi) + e + (p 2 ) -> /(gi) + /(gz) + n(j(kj)) + n'(j(kij). 

In the present description of the EEX matrix element and in the following sections on 
the Monte Carlo algorithm, we use an alternative notation for the fermion momenta: 

Pl=Pa, Ql=Pc ^ 

P2 = Pb, q2 = Pd- 

This is because such a notation has already been used for a long time in low-level MC 
programs and in the EEX distributions; we have therefore decided to keep it also in the 
relevant parts of the paper. We hope that it will help the user to understand the code 
more easily. On the other hand the notation with Pa,A = a,b,c,d, and ki,i = l,...,n, 
is very handy in building up a very compact tensor notation, with letters for fermion 
indices and numerals for photon indices, in the tensor notation for the multiple photon 
spin amplitudes; see ref. [jT2]| , the previous section, and the relevant parts of the program. 
Denoting the Lorentz invariant phase spaceQ, 

n j3 / 11 \ 

dr n {P; Pl ,p 2 , ...,p n ) = J] -Jj± *«> [P - 5>J , (14) 

the 0(a r ) EEX total cross-section 
°° 00 1 1 f 

-JZSZ^T^T / dr n+n'+2(Pi +p 2 ]qi,q2,k 1 ,...,k n ,k[,...,k' n ) py EX , r = 0,1,2,3, 

n=0 n'=0 ' ' 

(15) 

is expressed in terms of the fully differential multiphoton distribution 

p£L(Pi,P2, ft> ft> h, k n , k[, k' n ) = e Qlr*M^)+Q}^ F M 

n n' ( 

JI 25/(^0 Q(Q l] k,)]j2S F (kl)Q(Q F] k{)hj;\x,p 1 ,p 2 ,q 1 ,q 2 ) 

3=1 1=1 { 

+ ftij(^Pi»P2,gi,g2,fcj) + PiF ( X ' Pi' ^ 2 > ft' ft' k ti 



°EEX 



) Pi: P21 gi; g2; kj> kk) P^ff^-^i p^i gi; gz> ki, k n 



n> J>fc >i Srik^Siih) n ,> />m >! S F (ki)S F {k m ) 

| y^v^ P 2 t > F (X,p 1 ,p 2 ,q 1 ,g 2 ,k j ,k l ) | /^^(X, p 1; p 2 , gi, 92, fcj, fcfc, 

S I (k j )S F (ki) n >^>i>i ~Si{k 3 )~S I {k k )~S I {k l ) 



(16) 



2 This phase space has a slightly different normalization from that of the previous section. 
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The YFS soft factors for real photons emitted from the initial- and final-state fermions 
read 



. ^-<M«>' (17) 



The IR domains Qi and Qf are different for ISR and for FSR - this is easier to implement 
in the MC. We define f2/ with the condition k° < E min <C y/s in the rest frame of 
P = pi + p 2 called the CMS frame, and Qp with k'° < E' min <C 2m j in the rest frame of 
Q = qi + q 2 referrred to as QMS. The rest frame of X = P — Ylj=i n kj will be referred 
to as the XMS frame. In EEX we generally use small mass approximation m e ,mf, <C y^s 
(in the following sections we shall however use finite my wherever necessary). In this case 
the YFS form factors of eqs. flT2]) for the above JTs are very simple: 



V(Pl+P2) 2 4 71 \ 2 3 

2 a/ 2p ± p 2 \ o^2 a f, 2 ?i92 



7 = 7e = 2Q - Mn-^-1 , j f = 2Q}- hi^-l 

Again the IR-finite functions P^\k = 0,1,. ..,r (defined at the level of squared spin- 
summed spin amplitudes) are calculated, in principle, from the Feynman diagrams and 
they are defined completely in ref. \\12\f[ 

Here, let us only give the explicit definition of lowest order 

1 <7lT^ 0rn 

^\x,p u p 2 ,q u q 2 ) = - J2 ^-( X2 '^) ( 19 ) 

fc,/=l,2 

where 

0u = Z(pi, gl), t?i2 = Apu -Q2), $21 = A~P2, qi), $22 = A~P2, -fa), (20) 

and where all three-vectors are taken in the rest frame of the four-momentum X = 
Pi + P2 — Yuj=i called the XMS. The above /3q function is proportional to the Born 
differential cross section da Born (s,$)/dQ, for e + e~ — > ff. 

2.2 Spin effects 

Spin effects are treated in the most general way. For initial beams we use as an input the 
full spin polarization vectors and for the outgoing r we contract the r-production spin 



3 In reality they are constructed there by combining in a clever way first-order results and results of 
the triple convolution of the Altarelli-Parisi kernels. This solution has several practical advantages. 
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amplitudes with the polarimeter vectors of the two decaying r's defined in exactly the 
same rest frames of the r's where the spin quantization axes of the r's were defined in the 
first place. To find out exactly these frames was a little bit of an exercise, because for the 
calculation of the spin amplitudes ]ID| for the e~e + — > f fn{^) process we use the Kleiss- 



Stirling technique, for which these spin quantization frames have been found in ref. |TT 
We call them the GPS frames of the fermions (beam and t^). At the practical level, we 
have written a subprogram that performs the Lorentz transformations from the four GPS 
frames to the CMS frame. This subprogram is used to transform all r decay products 
from the GPS r rest frames to CMS frame. The polarization vectors of the incoming e 
have also to be provided in the GPS frames of e ± . It is simple to obtain by consecutive 
Lorentz transformation from the experimentally defined frames of e ± to CMS frame and 
later to GPS frames. It amounts to making the three-dimensional Wigner-Wick rotation 
of the beam polarization vectors from the experimentally defined frame to the GPS frame 
(in practice it is a rotation around the beams, due to the smallness of the electron mass). 
All of the above technique is an extension of the methodology used in KORALB |0|. For 
more details, see refs. [TO, III]]. 



Let us finally touch briefly upon another very serious problem relevant to implemen- 
tation of spin effects and its solution. In eq. @ the single spin amplitude 97tn ~* already 
contains 2 n (n + 1) terms (due to 2 n ISR-FSR partitions). The grand sum over spins in 
eq. (0) counts 2™4 4 4 4 = 2 n+16 terms! Altogether we expect up to N ~ n2 2n+16 operations 
in the CPU time expensive complex (16 bytes) arithmetics. Typically in e~e + — > fi~fi + 
the average photon multiplicity with k° > 1 MeV is about 3, corresponding to iV ~ 10 7 
terms. In a sample of 10 4 MC events there will be a couple of events with n = 10 
and iV = 10 12 terms, clearly something that would "choke" completely any modern, fast 
workstation. There are several simple tricks that help to ease the problem; for instance, 
objects such as £"c?t and the s-factors are evaluated only once and stored for multi- 
ple use. However, this is not sufficient. What really helps to substantially speed up the 
numerical calculation in the Monte Carlo program is the following trick of photon spin 
randomization. Instead of evaluating the sum over photon spins <X;, i — 1, n, in eq. (|3|), 
we generate randomly one spin sequence of (eri, o~ n ) per MC event and the MC weight 
is calculated only for this particular spin sequence! In this way we save one hefty 2 n factor 
in the calculation timej]. Mathematically this method is correct, i.e. the resulting cross 
section and all MC distributions will be the same as if we had used in the MC weight the 
original eq. (|3|) (see a formal proof of the above statement in Section 4 of ref. ||22|| ). Let 
us stress again that it is possible to apply this photon spin randomization trick because 
(a) the typical high-energy experiment is blind to photon spin polarization, so that we 
did not need to introduce the polarimeter vectors for photons in eq. @), and (b) for our 
(circular) choice of photon spin polarizations the cross section is rather weakly sensitive 
to them, so the method does not lead to any significant loss in the MC efficiency 



4 The other 2 n factor due to coherent summation over partitions cannot be eliminated, unless we give 
up on narrow resonances. 
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2.3 Electroweak corrections 

Electroweak (EW) corrections are implemented using the DIZET library ||23|| , more pre- 



cisely using its most recent version, 6.21, the same as is used in the ZFITTER [24 



semi-analytical program^. The method we use to implement EW corrections is rather 
simillar to that in KORALZ |l[: they enter through modified coupling constants in the 
Born cross section, and in order to save CPU time, correction factors are stored in the 
look-up tables; see the relevant section later. 

In the absence of EW corrections, coupling of two neutral bosons 7 and Z are defined 
in a conventional way: 

GZJ ZJ \ ZJ (~CiJ ZJ \ , _ D r 

A — 9v ~ a 9a 1 ^A — 9v 1 A — — K i L i 

glf = Q e = l, gf = Q f , gY = 0, gf = 0, 

z , e = 2T*-4Q e sm 2 6 w ZJ = 2Tj-4Q f sm 2 9 w (21) 
^ v 4 sin 9w cos 9w ' V 4 sin 9w cos 9w ' 
gZf = 2T? „ = 2Tf 

4 sin 9w cos 9\y ' 4 sin 9w cos 9\y ' 

where T 3 is the isospin of the left-handed component of the fermion (T| = — 1/2, T e 3 = 
-1/2). 

As in KORALZ [|l|, we implement the EW corrections as follows: the 7 and Z propaga- 
tors are multiplied by the corresponding two functions (scalar formfactors due to vacuum 
polarizations): 

1 C M 2 
H, = —— } H z = 16sm 2 9 w cos 2 9 w ^7= Pew- (22) 

In addition the vector couplings of Z get multiplied by extra form factors. First of all we 
replace 

z>e _ 2Tj-4Q e sm 2 9 w 2T e 3 - 4Q e sin 2 9 w F EW {s) 

9V - a - a „„„ n ^ 



4 sin 9w cos 6w 4 sin 9\y cos 9 W 

2Tf -AQ f sin 2 9 W 2Tf - 4Q f sin 2 9 w F EW (s) 



?f , , (23) 

9V - a - a „„„n ^ 



4 sin 9w cos 9w 4 sin 9w cos 9 



where F EW (s) and F EW (s) are electroweak form factors provided by the DIZET program 
and they correspond to electroweak vertex corrections. 

The electroweak box diagrams require a more complicated treatment. In the Born 
spin amplitudes we have essentially two products of the coupling constants: 

riZ.e^ZJ 1 Z,e \ Z,ew ZJ . x ZJ\ Z.e ZJ \ Z,e ZJ . \ Z,e ZJ Z,e ZJ 

G x G -x = \9v ~ X 9a )\9 v + X 9 A ) = 9 V 9 V ~ X 9 A 9 V + A 9 V 9 a ~ 9 a 9a > 

nZ^^ZJ 1 Z,e x Z,ew ZJ x ZJ\ Z,e ZJ x Z,e ZJ x Z,e ZJ . Z,e ZJ 

G x G x ' = \9v ~ X 9a )\9v - x 9a ) = 9v 9v X 9a 9v ~ X 9 V 9a + 9a 9a ■ 

(24) 



5 We would like to thank the authors of DIZET for help in implementing DIZET in the K.K. MC. 
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In the above, the following replacement is done for the doubly- vector component 

z, e zj 4T e 3 T| - STjQfFjwjs) - ST]Q e F% w (s) + 16Q f Q f F e J w (s,t) 

9v 9v == ' > 77 • 2 n 2o ' ^ > 

16 sin &w cos z Uw 

where the new form factor F^f w (s,t) corresponds to electroweak boxes and is angle de- 
pendent. 

In CEEX the Born spin amplitudes modified in the above way enter everywhere as 
building block, also in the spin amplitudes for the single and multiple real photons. In EEX 
the Born spin amplitudes with the modified couplings enter into the Born differential cross 
section, which is the basic building block in all EEX differential distributions. In both 
cases it thus constitute an improvement with respect to KORALZ, electroweak corrections 
are calculated for every occurrence of Born-like building block of the amplitude and not 
only once per event. Note also that we do not need to discuss final-state mass terms at 
this point. 

2.4 Beamstrahlung 

Beamstrahlung is the effect of loss of the beam energy due to the electron-bunch inter- 
action. It is essentially the emission of two effective photons strictly collinear with the 
beams, as an additional strictly collinear ISR. 

In the MC program, we require that the 2- dimensional beamstrahlung structure function 
be the "user function" , supplied or easily replaced by the user, without any loss of efficiency 
of the MC program. 

The most general form of the modification of the differential distributions due to 
beamstrahlung in e~(p a ) + e + (pb) — > X reads as follows: 



da bst (p a ,p b ; X) = J dz 1 dz 2 V(z 1 , z 2 , \/s)da(z 1 p a , z 2 p b ] X), 



(26) 



where T>(z\, z 2 ) is the double differential beamstrahlung "luminosity spectrum" for simul- 
taneous beamstrahlung from both beams. The luminosity spectrum is assumed to be in 
the most general form 

V(z h z 2 ) = 5(1 - Zl )5(l - z 2 )p + 5(1 - zi)pi(z 2 ) + 5(1 - z 2 )pi(z 1 ) + p 2 (z^ z 2 ), (27) 

where po is a constant and pi(z) and p 2 (zi, z%) are analytical functions for z e [0, 1). More 
precisely we allow for power-like integrable singularities at Z\ — 1: 

Pl (l- e )~e a , p 2 (l-e,z 2 ) -e^, p 2 (z u 1 - e) ~ . (28) 

The functions pi are regarded as completely arbitrary and we require that the /C/C MC be 
able to cope efficiently with any beamstrahlung luminosity spectrum of the above most 
general type. 

In the present version of the program, we include the beamstrahlung spectra as im- 



plemented in the CIRCE package of T. Ohl | 25| . The CIRCE package is based on the 
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results from the machine simulation program GuineaPig for linear colliders |26fl . Gener- 



ation of the Zi variables is done at the very beginning of the generation, together with 
the emission of the QED ISR total photon energy. The technical problem to be solved 
is hence the following: How to generate up to three variables, the two variables z\, z 2 for 
beamstrahlung and one z — 1 — v for ISR, according to an arbitrary, highly singular prob- 
ability density distribution? Not only are the beamstrahlung spectrum and ISR photon 
distribution strongly singular, but in addition the MC algorithm has to deal efficiently 
with singularities due to resonances and thresholds in the reduced CMS energy variable 

s" = SZiZ 2 Z. 

Furthermore, owing to the presence of the <5's in the general beamstrahlung distribution 



of eq. (|28|), four branches in the MC generation are present: one 1-dimensional, two 2- 



dimensional and one 3-dimensional, corresponding to each term in eq. (|28|) . In each 
branch the corresponding subset of the zi, z 2 and z variables is generated, using the special 
general-purpose MC generation program Foam |[27|| , developed recently for exactly this 
type of problem. In the present version, we also include another solution for the above task 



based on the classical Vegas program pq| , customized to our needs. This second solution 
based on Vegas is, however, much less efficient than the principal one (Vegas will probably 
be removed in the future version). Vegas is extremely efficient in calculating the value of 
an integral using the MC method, but our problem is different: the variables Zi, z 2 and z 
are generated at the beginning of the whole /C/C MC algorithm; consequently they have to 
be generated with the weight equal 1 (otherwise we would waste CPU time for calculating 
for instance the complicated QED matrix element for "unworthy" weighted events). To 
this task, Vegas is not very well suited (see a more detailed discussion in ref. p7j ). The 
new MC tool Foam is specialized for exactly this difficult task of generating unweighted 
events efficiently, according to an arbitrary multidimensional probability distribution. 

Even if we use powerful MC tool such as Foam, it is worthwhile to improve the effi- 
ciency of the MC generation by the appropriate change of integration variables (mapping). 
Let us take the most complicated case with J dz\dz 2 dz where z = s'/s, and s' is "after 
ISR". By the trial and error method, we have found that the best mapping for y/s < 1 
TeV is: 

xi = l-z 1 , x 2 = 1 - z 2 , v — 1 zx\ — r\ ,a , x 2 = r\ ,a , v = r\ h v m& ^ (29) 

where < r« < 1 are used by Foam, a is from CIRCE and, for the ISR, we have 
7 ~ 2(aQED/ir) \\i(s/m 2 e ). The integral gets transformed into: 

J dzidz 2 dzF(zi, z 2 , sziZ 2 z)Q{\ - Z\Z 2 z - v max ) 
f 1 VV X X 

= / dndr 2 dr 3 — — — -6(1 - z x z 2 z - v max )F(z 1} z 2 , sz x z 2 z) (30) 

Jo r 3 j r x OL t 2 ol 

f 1 , , , Wmax©(l - Z X Z 2 Z - V m3x ) 

= dndr 2 dr 3 - — ^j— — t F{z u z 2 ,sz 1 z 2 z) 
Jo [ax l ){ax 2 )(7^ 7_1 ) 

Note that this mapping only takes careQ of singularities at — 1; for the Z resonance 

6 Surprisingly, a more sophisticated mapping gives worse efficiency for both Foam and Vegas! 
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peak and the phase-space limit 1 — z\z 2 z < v max , the importance sampling is done by 
Foam (or Vegas). 

For energies up to 1 TeV the efficiency of the MC is better by a factor of more than 
10 with Foam than with Vegas. 
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3 Monte Carlo algorithm 



In this section we describe in detail the numerical Monte Carlo algorithm used to generate 
final-state four-momenta, i.e. points within the Lorentz invariant phase space, according 
to eq. (HI) for CEEX and eq. ([II]) for EEX. We shall not only describe the actual algorithm 
implemented in the present MC but also try to explain why we opted for certain solutions, 
and not for others. 

We start with a brief discussion of some important general issues. First of all, we 
would like to stress that we treat the MC technique as a means of integrating exactly 
over the phase space, without any approximationF]. Generally, our approach is that of 
the textbooks on quantum mechanics: the differential cross section is the phase-space 
times the scattering amplitudes from Feynman diagrams, squared, spin-summed. The 
MC technique is basically the numerical method of integration over the phase space. Our 
MC is in fact more than the phase-space integration, because we require that events 
(lists of four-momenta) be generated with weight equal 1, i.e. we directly simulate the 
scattering process. It is quite a strong restriction on the MC algorithm, see below, and 
it means that our MC program is not merely a phase-space integrator, but the full-scale 
MC event generator (MCEG). 

Concerning the Monte Carlo techniques, let us remind the reader that there are only a 
handful of elementary techniques such as weighting-rejection, mapping and multibranch- 



ingn - the whole art being to combine them into one bigger algorithm, see ref. [22] for 



more detailed discussion. We shall also generally follow the notation and terminology 



of ref. pUfl . There also exists a special group of self-adapting MC techniques/programs, 
such as the very popular MC integrator Vegas f28|| , which work, in principle, for arbitrary 
integrand distributions (just like a standard Gaussian integration program integrates any 
external function). We use in /C/C, for this role, the newly developed Foam [27| self- 
adapting program and the old Veskol (part of YFS2 ^0[) as building blocks of our MC, 
along with other elementary techniques. 

In general, we tried to minimize the use of the multibranching technique and rely as 
much as possible on the method of reweighting, constructing several layers of weights - 
the total weight being their product. Nevertheless, we have three multibranchings, one for 
the types of the final fermions / = e, /i, r, d, u, s, c, b, another one for the photon partitions 
(for n photons there are 2 n ways of attributing photons to ISR or FSR), and finally for 
helicities of the bremsstrahlung photons (for n photons there are 2 n helicity assignments). 
Each of these three multibranchings is well justified and unavoidable see below. 

In the development of the present MC algorithm there are certain critical issues, in 
other words, there are some problems that had to be solved, otherwise it would have 
been practically impossible to realize the entire MC. Typically, we have seen two possible 
solutions of an important problem, and we have opted for one of them. In the following 



7 It is necessary to stress it because in the MC event generators for QCD this may not be the case. 
In particular, the four-momentum conservation may be imposed at the end of the generation by certain 
ad hoc adjustments, which introduce an uncontrolled correction to phase-space normalization. 

8 The multibranching technique is referred to in other papers [E9[ as a multichannel technique. 
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we list these critical issues, explaining how the solution was found and/or the critical 
option chosen. 



Problem of the sum over partitions in the weight due to IFI. It was realized already 
in ref. |H] that the effect of the quantum-mechanical interference between photon 
emissions from initial- and final-state charges (IFI) can be introduced by means of 
the weight w IFl = p(q c ,qd,ki,...,k n )/p'(q c ,q d ,ki,...,k n ), where p includes IFI and 
p' neglects it. In MC events are generated according to p' in the first place. Both 
p and p' include the sum over 2 n ways of attributing photons to initial and final 
states, the so-called photon partitions. The natural way of generating p' is to use 
multi-branching and generate an event for a single partition at a time. Generation 
of p' involves the introduction of the additional "kinematical" weight of its own. 
The possible complication to be avoided is the sum over these kinematical weights 
over all 2 n partitions. Fortunately, this can be done: it is sufficient to multiply w lFl 
by the kinematical weight for a single partition, the one that is actually generated. 
The formal proof of this can be found in ref. |22] . The acceptance rate for n photons 



due to w 1F1 is 2~ n . It is a severe problem and only the partial, brute-force solution 
is applied at present. A better solution, taking into account that IFI is destructive 
for backward scattering and constructive for forward scattering, is needed in the 
future. 

Photon helicity generation. The sum over partitions in w 1FI already costs a lot of 
CPU time. Another sum over 2 n photon helicities would render the project imprac- 
tical, even with the present CPU processors. The solution is to introduce multi- 
branching for 2 n photon helicities and generate photon helicities event by event[]. 
It is justified, because high-energy experiment detectors are insensitive to photon 
polarizations; consequently the various photon helicity configurations (branches) 
contribute equally. Invoking the principles of ref. |22| , we can calculate weights like 



w 1F1 for just one helicity configuration, the actual one. The sum over 2 n photon 
helicities is avoided. 

Fermion-type generation versus photon- energy generation. In the MC we generate 
flavours according to cross sections that take ISR into account in an approximate 
wayp ! and the MC weight corrects it to the true value later. The approximate 
cross section involves a numerical integral over the ISR total energy, which has 
to be done before the MC generation starts. The integral can be done (i) for each 
flavour separately, (ii) or just once, for the flavour-summed cross section. In the MC 
generation the above corresponds to two options: (i) the fermion type is generated 
first, and the total energy loss due to ISR, just for this fermion, is generated next; 



9 Such a "randomization" of photon helicities was already used in ref. |2l| . 

10 Fermion-type (flavour) could be generated more simply, according to a crude cross-section equal to 
the Born cross-section, and the M C would correct the flavour rate in a way depending on the radiative 
correction, differently for each type of a fermion. Such a solution is unacceptable because, in the presence 
of the Z, ISR corrections are huge, a factor ~ 3, and the method would be rather inefficient. 
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(ii) the total energy loss due to ISR is generated first, for all fermions, and the type 
of the fermion using the Born cross section for the reduced centre-of-mass energy is 
generated next. We have chosen the second option. 

Comoving frame for the generation of FSR photons. Generating FSR photons 
is much more difficult than ISR photons because momenta of the final fermions 
"move", i.e. they are themselves variables in the phase-space integral, contrary to 
beams in ISR, which are fixed (except for beamstrahlung). Moreover, FSR emission 
distributions are also the simplest in the rest frame of the final fermions. For the 
construction of the FSR algorithm, it was critical to reparametrize the phase space 
in such a way that photon integrals are formulated in the reference frame attached 
to outgoing fermions. This comoving frame was used for FSR generation in the 
earliest version of the YFS3 MC which was incorporated in the KORALZ program 
since the version distributed at the time of the 1989 LEP1 workshop. The detailed 
derivation of the above phase space re-parametrization introducing the comoving 
frame can be found in ref. |32|| . 



Common IR boundary for ISR and FSR. Since ISR photons are generated in CMS, 
and FSR photons in the rest frame of the final fermions (QMS), the easiest is 
therefore to introduce the IR cut for real photons in terms of minimum energy 
in these two frames. This defines the IR boundary, and IR domains inside them, 
which are different for ISR and FSR real photons (the Lorentz-boosted sphere is 
ellipsoidal). This is not a problem, as long as the IFI is neglected. For the CEEX, 
however, the IFI is present and the IR boundary has to be common (for the basic 
MC in which IFI is neglected). The problem may look very difficult; however, an 
inspection of the proof of the independence of the total cross section on the IR 
boundary for the non-IFI case tells us that we may simply take the common IR 
domain, which contains both ISR and FSR domains, and, for each event, we may 
"remove from the record" all photons that are inside the new common IR domain. 
The above is true, however, only for events weight 1, not for the weighted events, 
that are present internally in the MC. For weighted events, such a procedure of 
"photon-removal" can still be done, provided it is accompanied by the additional 
weight, calculable analytically. The above procedure of the "photon- removal with 
the weight" was already implemented in the earliest version of YFS3 in KORALZ 
(although it was not really necessary for the EEX matrix element there); however, 
it was never documented^. 

Small photon mass as IR regulator in the MC. There exists all the time an interesting 
option in the MC implementation of the YFS exponentiation, which was never 
exploited - that is to use the photon mass as an IR regulator in the MC. The back-of- 
the-envelope estimate shows that its disadvantages are a slightly more complicated 
algorithm for the generation of the photon momenta and higher photon multiplicity. 



11 The analogous procedure was introduced in BHLUMI, where it was documented ||], and later on 
exploited for the construction of BHWIDE pi). 
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A clear advantage is that such a cut-off is Lorentz- invariant: this would therefore 
make it automatically the same for ISR and FSR, so that we would not need to do 
any gymnastics with the IR boundaries, as described in the previous point. 

The user who is familiar with our programs KORALZ and BHLUMI may be puzzled 
by the fact that we often we point out to some parts of the basic MC algorithm for QED 
bremsstrahlung, which have existed for a decade but are still unpublished or undocu- 
mented; let us comment briefly on this point. On the one hand, the full-scale ISR and 
FSR multiphoton generator YFS3 with the EEX matrix element was never published as 
a stand-alone program, for instance because of the lack of EW corrections, and it was 
absorbed into KORALZ. Its ISR algorithm was rather well documented in ref. [0O| , but 



its FSR algorithm was never in fact documented in detail. In the following sections we do 
this for the first time. Another reason for this time-lag is that the role of the MC programs 
at LEP1 for fermion pairs, perhaps with the exceptions of Bhabha scattering and r pair 
production, was limited to removing detector acceptance, and was not really a primary 
source of the SM predictions. The simple two-particle final state could be described fairly 
well with simple semi-analytical calculations such as ZFITTER. In LEP2 the role of the 
MC as a source of SM predictions seems to be more important, mainly because of the 
rise of the IFI corrections, and because of the very strong phenomenon of the Z radiative 
returnP]. 

In the following we shall describe the algorithm of the Monte Carlo generation of the 
events according to CEEX and/or EEX differential distributions. The algorithm is built 
using elementary methods of weighting and multibranching with occasional use of the 
mapping (change of integration variables) wherever it is possible and/or profitable. The 
weights are usually products of several component weights ordered in a "chain". Their 
role is to counter another "chain" of simplifications that were introduced in order to 
turn very difficult original differential distributions into simpler ones, in which we can 
integrate manually over certain integration variables (that is to generate them with the 
help of mapping), one by one. The remaining small subset of variables for which we are 
not able to perform the manual/analytical integration/mapping is treated with the help 
of the self-adapting Monte Carlo tool - in our case it is Foam and/or Veskol (optionally 
Vegas). The above "bottom-to-top" procedure of simplifications countered by weights, 
multibranchings, and integrations/mappings will be described in the following subsections, 
starting with removing IFI (countered by w 1Fl ), then simplifying the SM/QED matrix 
element, reorganizing phase space (mapping), integrating simplified ISR and FSR emission 
distributions (mapping), summing over photon multiplicities, etc., such that in the end 
we are left with only the integral over photon energy loss due to ISR and optionally due 
to beamstrahlung, which is fed into Veskol or Foam. In the MC program the order 
of action is "top-to-bottom", i.e. the photon energy loss due to ISR is generated first 
(together with beamstrahlung, if present), then the type of the fermion, the ISR photon 
multiplicity and momenta, the FSR photon multiplicity and momenta, and the series of 
weights brings us back to the desired EEX or CEEX differential distributions. The order 



12 It is even more true for the W W + channel because of the 4-body final-state kinematics. 
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of the multibranching and weighting elementary methods can be interchanged, with some 
care, and the rules for doing it are given in ref. p2| . We shall exploit this possibility, as 
already indicated. 



3.1 Weights and distributions in general 

Let us describe the whole organization of the weights and distributions in our MC. We 
have four principal distributions, pure phase space, model, crude and primary. Their 
ratios are the principal weights in the MC. 

• The pure Lorentz-invariant phase space (LIPS) distribution of eq. (0), which in- 
cludes four-momentum conservation 8's. It will not be generated directly in our 
MC0. It is our basic reference differential distribution in the sense that all other 
differential distributions of interest to us can be expressed in its terms: 

dcr(r u r n ) = p(n, r n ) dr n (P; n, r n ) (31) 

where the density distributions (we shall use this terminology) 

p(r 1 ,...,r n ) = da{r l ,...,r n )/dT n (P;r 1 ,...,r n ), (32) 

are analytical^ except for some simple discontinuities, with no 8's. 

• The model distribution is a density distribution corresponding to a physical model 

p Mod (ri, r 2 , . . . , r n ) = da Mod ( n ,r 2 , r n )/dr n+2 (P; t x , r 2 , . . . , r n ), (33) 

and it is our ultimate aim to generate MC events according to this differential 
distribution. Usually, we have to deal with several variants (perturbative orders) of 
the model distribution, with similar properties, that is to say, peaks. 

• The crude distribution is a density distribution 

p Cm (ri, r 2 , ...,k n ) = da Cru (r 1 ,r 2 , k u k n )/dr n (P; n, r 2 , h, k n ), (34) 

which is rather close to the physical model distribution, in fact it is close to all 
model distributions of a certain class, and it should be maximally simple. It should 
be Lorentz-invariant, with no unnecessary traces of the internal technicalities of 
the MC, just a maximally simple function of dot-products of the four-momenta. It 
should have the same pattern of peaks as all model distributions of a certain class. 
Here and later, r 1; r 2 will correspond to the momenta of the outgoing fermions, 
whereas ki will denote the momenta of the photons. In such cases we will write the 
phase-space dimension explicitly as n + 2. 



13 It is integrable, but not analytically, except for n = 2,3. 

14 We limit ourselves to the fermion pair production process. The r decays and quark-gluon parton 
showers are not included in the discussion. 
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• The primary distribution is a density distribution, 

dp P *%, £ n ), (35) 

defined primarily in the space S of variables £j with the following properties: 

— The integral f <ip Pn (£i, . . . , £ n ) is known independently; in most cases from 
analytical integration, or an independent numerical integration of the Gauss 
type, or from a special independent MC run. 

— A well-defined mapping r — > £ exists. The image of this mapping S LIPS does 
not necessarily cover the entire S. For £ e Slips the inverse maping £ — > r 
exists. 

We may therefore define 

p Pn ( r l> r 2, ...,kn) = rf<J Pn (£l, 6, • • • , in)/dT n+2 (P; T X , T 2 , ■ ■ ■ , fc„), (36) 

restricted to £ G Slips an d which is the distribution actually generated at the lowest 
level of the Monte Carlo. A weight will be assigned to MC points £ ^ Slips- The 
pPn correS p 0nc i s to events generated according to dp Pn , with all weights ignored (set 
to 1); it is "inelegant", with traces of many technicalities, not necessarily Lorentz- 
covariant, roughly similar to the physical model distribution. Its unique, great 
property is that it can actually be generated and integrated. Its integral f cr Pn sets 
the whole normalization of the MC. 

The main rationale for introducing the intermediate crude distribution, which obvi- 
ously stands between the primary and model distributions, is a very strong practical need 
of modularity of the MC program. For instance we want to use the same low-level MC 
event generator for both EEX and CEEX models. (As the example of BHLUMI shows, on 
top of the low-level with exponentiation one may even impose a model distribution with- 
out exponentiation.) We definitely want the MC event generator to have a well defined 
low- level MC part, which generates weighted events according to the crude distribution, 
that is with the weight 



Cm, 



W^(n, r 2 , . . . , r B ) = = £ E Slips, 

da FT1 (ri) p Fri {ri) (37) 

W Cru (r 1 ,r 2 ,...,r n )=0, ^ Slips, 

with the importance sampling corresponding to the entire group of physical models. The 
above weight is provided by the low-level MC to the outside world numerically, without 
any further details on how the event (ri,r 2 , . . . ,r n ) was actually produced - just as a 
black box. 

The model weight for the m-th model is the ratio 

T i/Mod/ ^ \ _ ^m" d ( r li r re) _ Pm° d ( r l, r n) /nn\ 

Wm (ri '-' rnJ " da^(n,...,r n ) ~ p^{ ri ,...,r n y ^ ] 
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calculated in a separate module, and in this module the crude distribution p Cru needs to 
be known functionally i.e., it is calculated locally, using a certain analytical expression in 
terms of the four- momenta of the event, and without any access to information from the 
lower-level MC (we even assume that such an information is already trashed). Of course, 
the total weight is the product of the two: 

W™ =W Clu W™ od , (39) 

and the total cross section is given by 

^ ot = <^r ) ^ Pri (40) 

This organization is definitely fully modular. 

Note also that although we did not require that the crude distribution be analytically 
integrable, such a property is strongly welcomed for the purpose of the technical tests. 



3.2 Crude differential distribution for EEX 

We define the crude differential distribution with respect to the standard Lorentz invariant 
phase space as follows: 

.Cru /„ „ . ;" I il ,J \ _ da Cru 



P[h,ri] (^1) Q2i ki, . . . , kn] k 1 , . . . , k n , 



dT n+n ' +2 (P; qi, <?2, fci, . . . , k h , k[, k' n ,) 



1 1 CTBorn(-Sx) S X 2 



Q yt j=l 1=1 

(41) 

where e e = 2E min / \/2p 1 p 2 , e f = 2E' mi J y /2q 1 q 2 /3 f = (1 - Am) / s Q ) 1/2 , s Q = 2q x q 2 + 2m 2 f , 
and we "dotted" the ISR photons in order to avoid a notation clash for CEEX in the 
next section. The above distribution features a maximum resemblance to the actual QED 
matrix element from the point of view of peaks and singularities, and is very simple. The 
infrared and collinear singularities are in S soft factors. The crBom^x) has a resonance 
peak in sx- The flux-like sq/sx factor is already present at 0{a 1 ) in the QED matrix 
element, and it can also be obtained from the leading-log approximation at any order. It 
is also necessary to include it in the primary distribution, in order to get reasonable MC 
weights at O{o 1 ) and higher orders. 

How did we get the above distribution? It was rather simple. We have taken the 
O(a ) version of eq. (|T5D, with only the [3^ term of eq. ([19]) in which we "averaged" 
over the angles in da Born and we have taken YFS form factors with the first term only, so 
that the IR finiteness is preserved. The kinematical factor sx/sq was adjusted to 0(a l ) 
and LL. The factor 2//3f is a pure convention. 

The above crude distribution is for EEX only; in CEEX it would only fit one single 
partition. 
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3.3 Crude for CEEX and multibranching over partitions 

In the EEX model we neglect IFI; we can therefore consider each bremsstrahlung photon 
to be attributed to the initial or final state; this is also true for the crude distribution 
for EEX. In CEEX the crude distribution of eq. (f|l]) is no longer usable, because it 
does not include the sum over partitions. To see it better, think about constructing 
the model weight as in eq. (|3q) , with the CEEX distribution in the numerator and the 
EEX crude distribution of eq. ([H]) in the denominator. The resulting weight would be 
wildly fluctuating. Let us therefore approach the problem in a more systematic way and 
construct a realistic crude distribution for CEEX by simplifying the density distribution 
in eq. (H), much as we did for EEX. Neglecting all spin effects (unpolarized case) and 
taking the O(a ) version of eq. @ we get 

p{p c ,Pd, h,...,k n ) 




rfLIPS„ +2 (P; Pc,Pd, 

k\, . . . , k n ) 

i „Y(C;p a ,...,p d ) P\(C>\ 1 

= _E I ZZl V" V" (Pklk 2 k n \ |W0) (pkikz fc n \l* 

7z' fluxfs) 4 ^ ^ \Aa"i(T2 ' ' ' cr n J \_ n Kx<n a 2 ■ ' • a„ J J ("42) 

V ' <r<=Tl A l=T l 

= ™h £ £ £ (rKi* <& *«) 5 1 ) (ikl* 

' M {p'} \i=l / \J'=1 

where s" = (p c + Pd) 2 - In the crude distribution we, of course, want to neglect the 
IFI. Neglecting it in the above formula means dropping non-diagonal terms p' 7^ p. In 
addition, in order to preserve IR cancellation we simplify the YFS form factor as well: 

Y(Q;p a , ...,p d ) -> 7ehi£ e + 7/lne/. 

As a result we obtain 

p(p c ,p d , A; 1 ,...,^) = - En^)E|*w| El®S^«)| a 7p&- 

{p} i=l (7, Aj 

(43) 

We can easily identify 

E |s^f| 2 = -8tt 3 ^(A; 1 ), S^ki) = Sjih), S (h) = S F (h), 

and also the Born-like expression 

X 2 Wn- Born 



A, 



which is weakly dependent on the dot-products s = 2p a p b , s" = 2p c p c , t = —2p a p c , 
t' = —2p b p d , u = —2p a p d , u' = —2p b p c ; it strongly depends on the X 2 in the Z resonance 
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propagator. We do not enter into the details of defining the above because we replace it 
by an "angular average" anyway: 



Aix 



Finally, the crude distribution for CEEX we define as follows 

p CrU (Pc,Pd, ki, k n ) = aa CEEX 

N C ' ' ' " " dr n+2 (P; Pc,Pd, ki,...,k n ) 



(44) 



1 1 rr Born ( X 2 \ X 2 9 n 

■ { P } Ht i=i 

(The overall factor like 2//3f, is just a convention.) It is almost obvious that the above 
distribution is the crude distribution for EEX of eq. (f|l]) summed over partitions - it 
is important to check it, at least to get statistical factors right. (We shall come to the 
question of the common IR boundary Q later on.) For example for n = 2, omitting final 
fermion momenta, we have 

Pf^ih, k 2 ) = I (2!pgfl(*i, h; ) + l!l!p£5(*i; h) + l!l!p^(A: 2 ; h) + 2!pg|( ; h, k 2 )) 

(45) 

and putting together [1, 0] and [0, 1] we have 

pf^(h, k 2 ) = pgjfo, k 2 ; ) + pf^ih; k 2 ) + pf ^( ; k h k 2 ). (46) 

As we now see, for arbitrary photon multiplicity we have the following relations between 
crude distributions for CEEX and EEX 

pffi(k 1 ,...,k n )= (47) 

n+n'=n 

where again we have "folded in" together all 

n\ n\ 



.n' / n\ n'\ 

contributions with the same ISR and/or FSR multiplicity. 

In the Monte Carlo we generate for a given n one of the distributions and we 

implicitly understand that in order to get all 2 n partitions we have to undo the above 
"folding in" , that is we have to perform the proper Bose symmetrization of the photon 
momenta; in other words, we must permute the momenta of all photons randomly, so 
that no trace is left of the primary ISR/FSR origin of a given photon. The above is 
always understood when we say that, we generate in the MC, all 2 n partitions using 
multibranching methods. 
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The above rule is important to remember, because the MC code itself may be mislead- 
ing for two reasons: first of all we exploit the equivalence principle of ref . |22| , in order to 



reorganize our weights, and the model weight for CEEX is multiplied by the crude weight 
for one of the partitions only, and not by the sum of crude weights over all partitions 
(branches) - so one may easily get the wrong impression that the multi-branching over 
partitions is absent. Secondly, in the low-level MC, when generating the primary primi- 
tive distribution, we do not generate first n and later n', with n = n — n' , but rather n' 
and h independently of each other, taking advantage of the particular properties of the 
Poisson distributions that govern them and allow us to do it. 

The bottom line is the following: the crude distributions for CEEX and EEX can be, 
and are provided by the same low-level Monte Carlo, because of the identity (f47|), and be- 
cause the multibranching over partitions is practically equivalent to Bose symmetrization. 
Understandably, this has great practical importance. 



3.4 Model weight and total weight 

The model weig ht for the 0(a^) EEX is the following 

WEEXWl) <?2, Kl, • • • , Kfi, fa,---, fa') - — 7Z-, • ; ~ 7T— , l 4J 



PeexIPuPz, gi, g2, fci-, fa, k[..., k' n , ) 
921 fa, ■ ■ ■ , k[, ... , k' n ,) 



where the model distribution in the numerator is from eq. flT6|), and the crude distribution 
in the denominator is from eq. (|4lD . 

The model weig ht for the 0{a^) CEEX is the following 

rrAr) , 11 1 \ P^C^EX.iPaiPbiPciPdi &lj ^2> • • • 1 hi) / Aa \ 

W^ EX {Pa,Pb,Pc,Pd, h, fez, ... , k n ) = — — ^ 3 ( n+2) _ 4 ; 49 

P[n] {PcPd, k 1 ,k 2} ... } k n )(2iry( n + 2 > 4 

where the model distribution in the numerator is from eq. and the crude distribution 
in the denominator is from eq. (H). The factor (27r) 3 ( n+2 )~ 4 is due to the difference in 



the normalization of dLIPS n and dr n . 

As explained in the previous section the corresponding total weight is 

WcEEx(Pa,Pb,Pc,Pd, h, k 2 ,..., k n ) = 

= W CEEx(Pa,Pb,Pc,Pd, h, k 2 ,-.., kn) W ClU (j) a , p b , p c , p d , fa, fa, ■ ■ ■ , k n) 



(50) 



and 



(51) 



Weex"^?!) 92; fa,..., k n ; k[, ... , k' n ,) = 

= Wbex(?i> 92; fa,..., k n ; k[,..., k' n ,) W° TU (qi, q 2 ] fa,..., k n ; k[, k' n ,), 

where W Cru is exactly the same, provided we do for CEEX the proper Bose-Einstein 
symmetrization. However, even this is strictly speaking unnecessary because the CEEX 
matrix element is Bose-symmetric anyway, by construction. 

Only one of many model weights is the principal weight used for a rejection of the 
events. Of course we choose the best one, that is the 0{a {2) ) CEEX-type. The other ones 
are available, and we use them for tests. 
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3.5 Phase-space reorganization 

Our starting point is the phase-space integral of eq. (HI]) for the crude total cross section, 
which can be rewritten as follows 



a 



oo oo 

Cru 



ds x J2J2 dr n+1 (P; k u ...,k n , X) -J[S e (kj)e e (kj) 

n=0 n'=0 J U ' j=l 

f dr nl+2 {X- ,k[,...,k' n , qi , q2 )±f[S f (kl)®M ^ j e^e^f 

(52) 

where P = p\ + p 2 . 

The integral above is Lorentz-invariant and, in principle, can be evaluated in any 
reference frame, not necessarily in the laboratory frame where P = andpi = (p°, 0, 0,p 3 ), 
which we call PMS. The implicit assumption is also that all momenta in the phase space 
are in the same reference frame at the time of the evaluation. The reality of the MC world 
is more complicated. Although the final product, the list of all four-momenta (MC event), 
is given to the user in one universal frame, the laboratory frame, in the intermediate stages 
the integral is split into Lorentz-invariant parts and each part is worked out separately, 
in the local Lorentz frame, which is the most convenient for generating subgroup of four- 
momenta in this subintegral. Then, momenta of this subgroup/subintegral are Lorentz- 
transformed to the laboratory. The MC generator is not fully documented if the relevant 
local frames and the Lorentz transformation connecting them to one another and to the 
laboratory system are not unambiguously defined. In this section we shall make an effort 
to do it. 

In the case of the above integral we take advantage of the Lorentz invariance of 
dr n ' +2 (X; , k[, . . . , k' n q%, q 2 ) and we transform all its variables to the frame where X = 
X = (yfsx, 0, 0, 0), the XMS frame, and put bars on top of them to mark this: 

a Cru 



/oo oo „ ^ n 

ds x ^2^2 dr n+i(P; h,---,K, X) — W_S e (kj)Q e (k 

n=0 n'=0 ^ U ' 3=1 

n' 

f dr n , +2 (X; , fcj, . . . , k' n q\, q 2 ) \ f[ ~S f {k[)Q } {k[ 
J n ' i=i 



, N O- Born(sx) Sx _2_ ^ z \ne e&V ln £/ 



4tt s q (3 } 

(53) 

This operation is still not defined unambiguously, unless we specify the direction of two 
space-like axes in the new frame or, equivalently, write down explicitly the Lorentz trans- 
formation Lx from XMS to CMS and back. 

Before we do this, let us note that there are at least three main ways of fixing the 
z-axis in the XMS frame. Two possibilities are to use as a direction for the z-axis three- 
momenta of beams, pi or — p 2 (they are not the same if s ^ sx)- We call them XMS1 



and XMS2 correspondingly. Note that these frames were introduced in refs. | 31| , |33 |. The 
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other choice is to use as z-axis in the XMS frame the direction of the P, that is the 
direction of the boost connecting PMS and XMS (that is the direction opposite to the 
total four-momentum of the ISR photons). This choice we shall call simply XMS, as it 
was used in ref. [13 j. A third choice is to do a so-called parallel boost Bx along the 



direction of the X in PMS. The corresponding transformation matrix is 



B 



x 



x° 



X T 
M x 

X®X 



X 2 



M 2 X , 



(54) 



M X (M X +X°)_ 

where T marks the matrix transposition and ® marks the tensor product. This is our 
choice in the present version of the program (the same as in the YFS3 MC). Note that, 
in general, in the XMS frame defined with the above transformation, the z-axis is not 
parallel to p\ or — p 2 . The transformation from our XMS to CMS is 



k' 



CMS 



Qi\CMS 



^xqi 



L 



x 



B 



x- 



(55) 



The transformations for other types of XMS are given explicitly in ref. |32 . 

Now comes an important point. As already indicated, the emission of the FSR photons 
is done in the comoving frame attached to the momenta qi of outgoing fermions, that is 
in the frame where Q = qi + q 2 = and q\ = (g°, 0, 0, \qf |). We shall call it QMS. So why 
do we not transform immediately from XMS to QMS, and generate photons there? The 
problem is that Q = X — K> an< ^ ^ n or der to get from XMS to QMS we have to know 
k\ in the first place. We are stuck. The solution is to reparametrize the FSR integral with 
the help of the integration over the Lorentz group; the details are given in ref. ||32|| . Here 
we just apply the result of this work and obtain the new formula: 



a 



Cru 



oo H n 



n=0 

dip d cos ui 



d 3 k, „ ~ 



„, jlx ./i-'S, (/,• )(-),(/,• ) 5[s x -(P -$>) 2 1 ' - ,u 



n=0 j=l 1 
C"Born(-Sx 



3=0 



E ^ / ds Q II ^ 25/(^)6,® s(s x -(Q + Y: k'j) 2 e^'t, 

n'=0 ' J 1=1 ^k i \ j =0 J 



(56) 

in which those variables with a tilde are defined in the QMS. Note that the Jacobian due 
to the reparametrization of the FSR integral cancels exactly the factor (sx/ sq)(2/ flf). 
The explicit transformation from QMS to XMS defines the meaning of the new ip,u 
integration variables: 



ki — L^ki, qi — Ly^qi, 

L A = R 3 (iP)R 2 (u)Bz\ X = Q-Y J ~k'y 



(57) 
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The most important fact is that the explicit integration over q\ and q2 has disappeared 
completely! The angles ip, to parametrize the three-dimensional orientation of the mo- 
menta set (k[, . . . ,k' n , qi, q 2 ) as a whole in XMS. 

The above treatment of the FSR phase space is the simplest one - this is why we 
adopted it. Note that the angles tp,u have no direct geometric meaning of polar angles of 
a certain momentum in a certain frame - they are just parameters in the Lorentz trans- 
formation. In particular, momenta q~ x and —q~ x do not coincide with the z-axis in XMS. 
Such an arrangement is, however, possible. It corresponds to a different transformation 
La and the integral fl56|) would look slightly different. Such an alternative solution is 



described in detail in ref. 32 



The crude integral of eqs. (f|l]) and fl56|) can be rewritten as follows 



oo oo 

Cm 



Y YY / dT n+nl+2 (P;qi,q 2 ,ki,...,k n ,k' 1 , 

t A ti a r> h n =C) « f — fl 



a 

f=/j,,T,d,u,s,c,b n=0 n'=0 

„Cru 



k' n >) 



Y / dSx a Lm( S x) / 



x P[nji']{q^i q^'i ki, ■ ■ ■ , k n ; k 1} . . . , k n i / 
dipd cos uj 



f=H,T,d,u,s,c,b " (58) 



oo ^ n 

77 n 



rl 3 k- ~ 

3 5 e (%)9 e (A ;i ) 5 I s x - (P 



n! ll jfcO »eK~,W 3 
n=0 j=l J 



oo ^ 
n'=0 



J r _i k\ \ 




This clearly factorizes into independent ISR and FSR parts, with the integration over 
the effective mass sx pulled out as a principal integration variable. The above integral is 
ready for the MC generation. 

Let us stress again that in all the above reorganization of the phase space for the 
cr Cm , from eq. fl4lf) through eq. (|52D to eq. fl58D , we only changed variables with no 
approximations, with full control of the Jacobians of all mappings. For completeness we 
write the total transformation from QMS down to CMS: 

k'i\cMS — BxL^k'^ qi\cMS — Px^Aq.%-, 

L A = R z {iP)R 2 {u)B~ jc \ X = X\ CMS , X = Q~Y~k', (59) 

j 

and X is defined in the rest frame of the outgoing fermions. 

In the following sections we shall introduce variables that are used directly in the 
MC generation, separately for ISR and FSR, and we shall define primary differential 
distributions generated at the lowest level of the MC algorithm. We start with the FSR 
- in this case the corresponding primary differential distribution sums up to unity, and 
next we elaborate on the case of ISR. The case of FSR will be described in detail, while 
the case of ISR will only be summarized, since it was already discussed in ref. 
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3.6 FSR momenta 

In the following we shall describe the MC algorithm for the generation of the FSR photon 
momenta. Although it was used for almost a decade |I|,|33|, it is the first time that it is 
described in full detail. Let us consider the FSR part of the crude integral of eq. (|58l): 

i B f n ' r d 3 k' ~ / / ~ \ ; 

',9 i=i ^'7 V ^ 1=0 ' 



4m > 



P f = (l-^ 2 , aiJ = — • £ /=^> Q = (v^, 0,0,0), 



where we restored finite mj, photon momenta k[ are defined in the QMS rest frame of 
the outgoing fermions - the natural reference frame to describe the emission of the FSR 
photons - and -E^in ls ^he real photon minimum energy in this frame. Let us now express 
the photon momenta in units of as well as introduce polar parametrization and 

other auxiliary notation: 



k\ = kj = Xj (1, sin 9j cos 0,-, sin 9j sin </>,-, cos 0A, 
J 2 2 



Z=0 

With the help of the above we are able to eliminate the 5-function: 



(61) 



/ 



Q(s Q (k u ...,k n >) - 4m 2 f ) 



Ava 2 j x / im^ 



(62) 



1 + K° + \Ri 



and from now on 



S Q = SgCh, M = 1+ / X + 1 R 2 - ( 63 ) 



Also the single-photon distribution gets transformed 

flu » n „. / ^,2 



{ m 2 \ 1 + B 2 u 2 1 u 2 1 
■H ) = ~~x~x ~ n I2 — o IT' % = l-/3/cos^, 8 ij = \+p J <xxO j 



s Q J 5^ 2 2 5 2 V 

(64) 
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and the whole integral is transformed into the semi-factorized form: 



OO 2-7T 1 



^ej -1 

The reader should not, however, be misled by the apparent simplicity of the above integral 
- it does not factorize yet into a product of independent integrals, one per photon, because 
the collective dependence on all photon momenta kj is entering everywhere through the 
variable sq, see eq. (|63|). 

Another complication due to the use of f */sq as an energy scale is that, in the case 
of the hard FSR photon, the upper limit of Xj extends to large values, not really to 
infinity because of the 0(sq — 4m 2 ); nevertheless, this is not very convenient for the MC 
integration. This is cured with the following change of variables: 



1 • ' 1 - E Vj ' 

1+Elj= 1 =1 + ^ = 1 + B^ = f3Lf 1 _^ (66) 

Y 1 - 22 j Vi S Q S Q \ s x 

which leads to 

n'\ 7 Vj J 2ix J 7T V s Q , 

o -l (67) 

1 + 7?° 



1 + K° + ±K 2 



6(sq -4m 2 ,) e^ ln ^ 



4 

2 



With the new variables the condition sq > Amj (easily implementable in the MC) trans- 
lates approximately into 2~2jVj < !■ Furthermore, we have 



1+& <1 



1 + K" + i/? 2 



which is ideal for the MC 0. The new IR limit y 3 - > £//(l + A' ) is however inconvenient 
for the MC. The solution is to substitute 

e, = 5 f (1 + K°) (68) 

where 5/ <C 1 is from now on the new IR regulator for FSR real photons. Note that this 
sets 



15 In the case of ISR a similar factor is causing a lot of trouble because it has a negative coefficient in 
front of K 2 and the corresponding contribution to the MC weight is not well bounded from above. 
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as a lower limit for the photon energy in the QMS, which is higher than the previous one 
-^min = \y/sQ^f (for e = 8f). Consequently, we have to keep the value of 5f very low, in 
fact we needj^] 8f <C mj-fsx, which can be a problem for / = e. 



3.6.1 Simplifications and MC generation 



Up to this point, the FSR integral of eq. floPf ) was transformed without any approximations 
and the integral was conveniently parametrized for the MC generation: 



2n I 



n ' 1 fJ[J Vj J 2n J 7T \ Sq J 

J " 1 <5 / o -l (70) 



' - + Q(s Q - Am}) e 7/M*/(i+^)). 



1 + K° + 



There is also a one-to-one correspondence between the points in the Lorentz-invariant 
phase space and the points in space of our new variables: 

{n',(k' 1} ... } k' n ,)} «-> {nf, (yj,9j,<f>j), j = l,...,n'}. (71) 

We can also write explicitly the differential distributions in the two equivalent parametriza- 
tions 

ddn' _ ®(s Q -4m 2 f ) 7/ ln 



r[25 / (fc;.)0(^-^min), 

ds Q 5(s x -{Q + Y.K))Yl^ 

^ 1 = ' j = l AK 3 

d$n> _ Q(SQ-Am}) elfHSf{1+R0)) ( a y* e( yj -6 f ) ( 

dy j d cob 0jd4>j ^ V 271-2 / f = \ v, \ j, sq)' 

(72) 

We are now ready to introduce the simplifications leading us to a primary distribu- 
tion, which can be integrated analytically and generated using standard uniform random 
numbers. The simplifications are[] 



/ ml\ _/ m 2 f \ 



i+pj i 



P f l-pj cos 2 Oj 

=r- - Q(sn -Ami) -»• 1, v 7 

1 + ^° + ±K 2 V Q V 



^flnidfil+K )) e 7/ln(5/) 



16 This should be listed as a disadvantage of the actual method of the MC treatment of the FSR. 

17 We drop the mass term from f(0j) for the same reasons as in the case of ISR; see next subsection. 
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where 

& f = (1 - {m)/s x Tyl\ 7/ = Q) I ^ In \±h . (74) 

The main purpose of the above is to remove any complicated dependence on the momenta 
of all photons through sq - it is achieved trivially by replacing sq by sx- With this hard 
FSR photons, get stronger collinear peaks at cosOj = ±1 in the primary differential 
distribution. The resulting FSR primary differential distribution is: 



IYj=i d Vj d cos 9j 



and the compensating weight transforming the primary distribution into the crude distri- 
bution is 



fsr ~ <r ~ i + r° + \k* u[sq 4m ^ j e 11 



W 

//.«?'•'•' 1 l /, " ..- 1 A-i • 



(76) 



Events {n', (yj, cos 9 j,<fij),j = 1, ...,n'} generated according to d^r m , defined in eq. ( [77] ) 
below with the weight u>pg R , will be distributed (if rejection is applied) according to the 
differential distribution, that is the integrand in eq. (ff2|), as desired. There is only one 
thing to be remembered: removing Q(sq — Am'j) in the primary distribution means that we 
cannot map every event generated according to d^J 1 into a Lorentz-invariant phase-space 
point, the unambiguous mapping exists in a strict sense in only one direction: 

{n',(k[,...,k' n ,)} -> {n'^yj^j^^J = l,...,n'}. 

This is, however, not really a serious problem because it occurs only for the events with 
w fsr = 0' while for events with w^sr ^ we are able to map 

{nf, (k[, . . . ,k' n ,)} <- {n',(yj,9j,<l>j),j = 1, . . . X}. 

Nevertheless, one should keep in mind that some 0-weighted events generated by the FSR 
algorithm do not have four-momenta assigned to them. 

Finally we may check that, as advertised, the integral over the FSR primary distribu- 
tion can be evaluated analytically: 



2 * 1 ,2- 



_ 5 / - 1 (77) 



oo 



n'=0 - J / n'=0 
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The photon multiplicity for the primary distribution is the standard Poisson distribution, 
with the average 

(n')=7/m^-, (78) 

and the overall normalization is trivially equal to 1, which is a natural choice for the FSR 
anyway. 

The MC generation of the distribution ( [To]) is rather easy. It is fully factorized - 
variables cos9j, (fij and yj can be generated independently. The distribution of (fij is just 
flat and the distribution of yj is trivial to generate: 

J = 2vrr 1 „ yj = S r / j , (79) 

where rjj are the standard uniform random numbers < < 1. The distribution of 
cos 8j requires applying the branching method: it is split into two components 

2 1 1 



l-^/cos 2 ^ 1 -p f cos 9j ' l + /3 /C os0/ ^ 

and cosOj is generated according to one component, chosen with the equal probability 
between the two. For example, if the first component 1/(1 — j3f cos 9 'j) is chosen then 

^H^-^^rrfr}' (81) 

where r^j is another uniform random number. 
3.7 ISR momenta 

In the following we shall describe the MC algorithm of the generation of the ISR photon 
momenta. The algorithm was already describedF] in ref. []30|], and for the sake of com- 



pleteness we shall describe it here, but without going into the fine details. Let us consider 
the ISR part of the crude integral of eq. (EBi) for one final fermion type /: 



I g7e In E e 



(82) 

where E m[n = e e \\fs is the minimum energy of the real ISR photon in the laboratory 
CMS. In the first step we introduce the variable v = 1 — sx/s and order energies of the 
photons 

{ ^ J j (83) 

6(fc? - k° 2 ) Q(k° 2 - A; 3 °) . . . Q(k° n - E min ) s(v - 2KP ~ ^ 



IS 



The essential part of the algorithm was given in ref. 36 
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where K = Y^j=o and t> max = 1 — Amj/s. Now we rescale all momenta and introduce a 
polar parametrization 

ki — rjki = ?7Xj(l, sin 6*4 sin0,, sin^j cos0j, cos^j); (84) 

we fix the scaling factor r\ such that k\ — x\ — v. 



- k° 2 ) <S>(k° 2 - k° 3 ) ...e(k° n - E mhl ) 5[v- — + — ) e"* I]l 



X=J drjSiv-ki/v) J dv al om (s(l - v)) f[ J ^S e (k,) 

KP K 2 

~ + lTj~ (85) 

"max Ti l 2 1 

= J *>oL.M-v)) II/ fj f 
o J -1 
<%-xi) 9(xi -z 2 ) 0(x 2 -x 3 ) ...0(A o x n -e) e 7eln£e J{K,v), 

where r/ is the solution^] of the equation v — "^j-?] + ^j-i] 2 = and J{K, v) is an overall 
Jacobian factor: 



^{K,v) = ^- - W - * =1(1 + 



(86) 



v/i f 2 v/i _ r ? 2 P 2 K 2 

Vo = ^k~o i + v /r^^ Ao ' A= (^ = (^- 1 ' ^^^^ 

and the photon angular distribution is governed by 

_ 2 2m 2 1 2m 2 1 

^ COS j > ~ (l-/3cos^)(l + /3cos^) ~ s (l-/3cos^) 2 ~ s (l + /5cos%) 2 ' ( ^ 

3.7.1 Simplifications and MC generation 

Up to this point the ISR integral of eq. (]82l ) was transformed without any approximation 
and we maintain, modulo ordering of the photon energies, the one-to-one correspondence 
of the points in the Lorentz-invariant phase space and in the space of the new variables: 

{n,(ki,...,k n )} ^ {n, (a^flj,^), j = l,...,n}. (88) 

Before we define the primary differential distribution to be generated in the MC, let us 
write once again explicitly the two equivalent (modulo energy ordering) parametrizations 



Note that for a single photon A = and rjQ = s 1/,2 /2. 
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of the ISR crude differential distribution: 



2fc° 



d3 n j ( /1 ^ f a \ s/ N @(Ao^ n — e) 



(if nj=i cos 0jc£0j Born \2tx 2 J x n (gg) 



n-1 



The simplifications leading to the ISR primary differential distribution are the following: 

2 



f (cos 6j) -> /(cos 6^) 



[1 -p cos 9^(1 + (3cos9 j ) 



^•)-*«-K i+ ^> (90) 

6(A x n - e) — > e(x n - e), 



where 



7 e = 2(a/vr)ln(s/m2). (91) 
The resulting ISR primary differential distribution is 

" 3 " p " -«i„Mi-«))C£Y«<«-*.) e(x "- e) 



rii=i dxjd cos 9 jd(J)j orn V^ 71 " 2 / 

n-l „ / _ x n 

n j g3 ~" ij n/( cos ^) e7einee n > °> (92) 



3=1 " J 3=1 



'0 



^Bornt 5 ) <*(*>)> " = °- 



and the corresponding weight is 



...Cm _ d ?n _ r\f\ _ f) J(K,V) 1-r /(cosflj) 

WlSR " aW ~ ( ° n } Jo(v) iJ/(oos^)- (93) 



Let us explain and justify the simplifications of eqs. (|90|). The replacement f(cos9) — ► 
f(cos9) is not really necessary in the present context of building an efficient MC algorithm 
for ISR. We could do without it, because the f(cos9) distribution is rather simple. The 
problem is really in the model weight at 0(a^ and higher orders. As is well known in the 
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soft limit, the helicity-non-conserving spin- amplitude contribution vanish, on the other 
hand, a perfect helicity conservation contradicts the angular-momentum conservation for 
a photon emitted exactly parallel to the respective emitting fermion, and this is reflected 
in the photon distribution 

/(COS 0y) = [(l_ /3cos ^.)(l + /3cos ^.)]2' ( 94 ) 

which has explicit zeros at cosOj = ±1. When we admit the exact Ola 1 ) hard photon 
emission matrix element, then for hard photons these regions close to cosOj = ±1 will 
be filled in by the helicity-non-conserving contributions, and the model weight based on 
/(cos#) would fluctuate wildly when we approach in division by /(cos#), division by 0. 
The solution is not to have these zeros at all, at the level of the primary distribution, 
and this is why we opted for /(cos#). In such a case the product of the model and crude 
weight will be regular at 0(a l ) and beyond. 

The other two simplifications are introduced for purely technical reasons. The simpli- 
fication J — > Jq is especially costly in terms of the MC efficiency because, for v — > 1, it 
introduces the "spurious" singularity (1 — v ) _1 / 2 . Together with the usual (1 — v ) _1 from 
<J Born( s (i — v )), it builds up strong singularity (1 — v)~ 3 / 2 in the primary differential dis- 
tribution, and huge primary integrated cross section J 1 4m // s dv(l /rrif. It 
is almost completely compensated by the very high rejection rate of events close to v — 1 
due to the ratio J^(K, v )/ J$(v) in The rejection rate is ~ (rtif/s 1 / 2 ) \n(Am 2 /s). For 
muons 99% of events are rejected. However, in most of the applications, this effect can 
be easily eliminated by setting v max = 0.999 or lower. 

The above problem is unfortunately unavoidable in the actual MC algorithm for the 
ISR. It can be traced back to the fact that the present MC algorithm is not very well 
suited for the emission of the two photons of large effective mass, such as the simultaneous 
emission of two hard photons along two beams. The present algorithm "folds in" together 
the energies of photons emitted from both beams. A more sophisticated algorithm, in 
which photons emitted from two beams are generated independently, is needed in order 
to eliminate this problem. 

The 9(Aox n — e) contribution to the weight has interesting consequences. As discussed 
in refs. [|30],|3(J it leads directlyQto a characteristic factor -F(7 e ) = e~ Cle /T(l + j e ) in the 
da/dv and in the total cross section. 

20 As shown in refs. |50| this is not completely straightforward. In fact F(j e ) is not present for v < e; 
it is nevertheless present in the integrated cross section, since the corresponding negative contribution is 
located just above v = e. 
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Finally, let us integrate analytically the ISR primary differential distribution 

„_n n 7 „•_! 70 X j JO Z7F 7-1 " 



J 

n=0 n=0 q j=l 

o7e Ine 



o x i 

5{v - x x ) e(x! - x 2 ) e(x 2 - x 3 ) . . . 0(x n - e) ln£e Jo(^) 



"max / qq 

J dv aLnW ~ v)) Jo (v) e> ln - Uv) + 6(« - e)-^ T^T)! (^ ln f) 

"max 



n-V 



(95) 

How do we generate the primary differential distribution o?J Pn ? We start with the 
generation of v according to 

-fo- = 4orn(s(l ~ V)) J (V) . (96) 

This is done by using a general-purpose MC tool such as Veskol or Foam; care is taken 
of any possible resonance or threshold in the (J^ orn (s(l — v)). In the next step, photon 
multiplicity n is generated. For v < e we have simply n = 0, and for v > e the photon 
multiplicity distribution is: 

?l" = const x -__U e ln-J , (97) 

which is just the shifted-by-one Poisson distribution P n _i, with the average Z.n — 1) = 
7 e ln(f/e). The angles cosdj and <j)j are generated in the same way as in the previously 
discussed case of FSR. 



3.8 Getting common IR boundary for FSR and ISR 

Let us consider the case of EEX: 



'eex 



W£2x da 



Cru 



(9* 



where the model weight W^x * s defined in eq. ( ^8|) in terms of the 0(a r ) EEX differential 
distribution of eq. ([TjjJ). Using eq. (|5^) with the later substitution 



,• | | , 2QK'\ 



(99) 



i=0 
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which was introduced in order to facilitate the MC generation, we obtain 

OO CO « 

a EEx{ A } = ^Y1 / d(7 £n'](^>^) A(n,k 1 ,...,k n ;n',k' 1 ,...,k' n ,;p i ,q i ) 

n=0 n'=0 

x ^ / EEx( n ' k±, . . . , k n ; n , k±, . . . , k n ,;pi, qi) 

dafZM,Q F ) = ds x ffB ° r ; W rfr„ +1 (P; A*, . . . , k n , X) e 7eln£e — TT 2S e (k j )G(tt I , kj) 
5 47r n! - LJ - 

i=i 

x dr„, +2 (X; fci, . . . , k' ,, q u q 2 ) ± e 7 ' 1 " -L TT 2S f (k{)Q(Q F , k[) , 

S Q Pf n l t = \ 

(100) 

where for the sake of the discussion of the IR cancellations we have introduced a general 
acceptance function A. Every physical, i.e. IR-safe, observable corresponds uniquely to 
one or more such acceptance functions. Just to give an example: the total cross section 
corresponds to A = 1, the forward-backward asymmetry is related to A expressed in 
terms of final fermion momenta like A = Q(qf), the cross section for the production of 
exactly two photons above E — 1 GeV corresponds to A = £\ . Q(k® — E )Q(k < j — E ), 
and so on. The acceptance function corresponding to a physically meaningful, IR-safe, 
observable has to obey one important rule 

lim A(n, k u . . . , fci-i, k u k i+1 , ...,k n ) = A(n - 1, k u . . . , k^, k i+1 , k n ), (101) 

and a similar rule should hold for FSR photons (in EEX we can make a distinction between 
ISR and FSR photons because we neglect ISR-FSR interference) . 

So far we kept the IR domains different for ISR and FSR, for ISR f2j was defined by: 
kj < Selv^ m the laboratory CMS system where p\ + f>2 = 0, while for FSR tip was 
defined by fc'J < S f ((s Q + 2K'Q) / s Q )\^fs^ in the QMS system where q x + q 2 = 0. Our 
task is now to bring the two IR domains together. 

We know that the total cross section and any IR-safe observable are completely inde- 
pendent of Q F and The self-suggesting solution is, loosely speaking, to set 5/ so small 
that we always have £l F C and simply neglect all FSR photons k[ e 50, = Oj\O f , 
that is just remove them from the list of the generated momenta in the MC. Note that 
because (sq + 2K'Q)/sq ~ sx/sq <C sx/i^mj) we may need Sf/e <C s/(4mj). 

Let us work out the details of the above method, providing a formal proof of its 
validity. The above prescription definitely leads to a certain new crude distribution 

in which the IR-domain O = Oj is common for ISR and FSR photons. The question is: 
What is the above new crude distribution? 

It turns out to be calculable analytically. (In the calculation we follow closely the 
algebra of the formal proof of the independence of the physical observables on the IR 
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domain Q as given in ref. ||12|| .) Let us consider the internal FSR subintegral in eq. ( |10C| ). 
that is all ISR photon momenta are fixed; 

oo „ 1 n' 

= W d ^+2(X; k[,..., k' n „ q u q 2 ) — J] 2^(^)6(0^, k[) b(k[, . . . , k' n ,; Pi , Qi ) 

n'=0 ^ ' Z=l 

X ^ / EEx( ri ' ^1; • • • 5 ^ni n ) ^1? ■ ■ ■ > ^n'SPi) (Zi) &1> • • • 5 ^ni n , k l7 . . . , k n ,]Pi, q^j. 

(102) 

Following f2/ = fi^ (J 5fi we split every photon integral into two parts and reorganize the 
sum factorizing out the integral over Sfl: 

^ = E ^ II { / -pr © W *0 + / y e(n„ fc|) SfW) 

n'=0 1=1 ' ' 

y dr n , +2 {X, k[] q x , q 2 ) b(k[, k' n ,; Pi , ft) 

e^e(:'){/§^^^)}- 

n'=0 s =0 V 7 W ^ 



/s n — s 

^„/ +2 _ s (X-^;fc / 1 ,...,/c^_ s ,g 1 ,g 2 ) JJ 0(ftj, £/(*;,') &(&i> • • • > *4'-*;P»> ft] 
1 i=i 

(103) 

where 0(5f2, fc') = 1 for k' G <5f2 and = otherwise. The most important ingredient in the 
above algebraic transformation was that the model weight W^ x , due to the particular 
expansion of p^ x into /3-components, see eq. (|IED, also fulfils the "IR-safeness" condition 

lim W E ft X (n', k[, . . . , k\_ Xl k' i: k' i+1 , . . . , k' n ,) = W^ x (n' — 1, k[, . . . , k[_ l} k' i+1 , . . . , k' n ,), 

(104) 

and consequently the function b(k[, . . . ,k' n ,;pi, qi) as well. The resulting integral 

oo „ 1 n 1 

W=E / dr n/+2 (X ] k[,...,k' nl ,q 1 ,q 2 )—]j2S f (k' l )Q(Q I ,k' l ) 

n'=0 J U -l=l (105) 

exp ( / w &m k) 2 ^ f{k) ) " " k ' n ' ] Pl} qi) 

gets an additional exponential factor, which is easy to interpret. It can be expressed in 
terms of the function 
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as follows 



6XP (/ l^ 6 ^'^ 2>§/(A;) ) =GXP ( 2 Q 2 f a B{n I ,q 1 ,q 2 )-2Q 2 f aB{n F ,q 1 ,q 2 )). 

We have therefore found out by explicit calculation that in the proposed method, in 
which for Qp much smaller than Qj we skip photons that fall into Sfl = Qi \Q F , the 
distribution of the remaining photons is the following 

dtf n %(n I ,{} I )=ds x aB °^ Sx) dr n+1 (P; h,...,k n ,X) e^^-^U^M^i^) 

71 U ' i=i 

2 1 n ' 

x dr nl+2 (X; k[,..., k' n „ q u q 2 ) -± - e R ^ — \[ 2S f {k[)Q{^ k[) 

s Q (3 f n'\ 

Rf = 7/ In fs f SQ+ s 2K ' Q ^ + 2Q)uB(tt h q u q 2 ) - 2Q 2 f aB(Q F , q 1} q 2 ). 



(106) 



Also note that, by construction, the integral value is preserved 



n,n' n,n' 

Now, since the IR boundary in the above new distribution has changed for FSR pho- 
tons, we cannot continue to use the pQ EX of eq. flTBl) . We have to use another p* FFX i n 
which we replace B(Q F ) by B(Qj) in the YFS form factor, see eq. ([TT]) : 

*( r ) _ ( r ) 2Q?a(B(n /) g 1 ,ga)-B(n F) gi,g2)) 
r EEX ~~ r EEX ' 

and consequently, since the model weight is the ratio of the model and crude distributions, 
the new exponential factors cancel out, and the new model weight is functionally exactly 
the same 

VV EEX — VV EEX' 

In the new MC calculation 

4L = j W£& d* c ™* = f Wj&x da c ™* (107) 

both the product of the weights and the normalization is the same. The above result is so 
trivial that, in fact, in the MC program for the EEX model we change almost nothing - we 
are only omitting hidden photons in the calculation of the model weight. This simplicity 
reflects the basic fact that very soft photons are unimportant for all IR-safe integrand 
functions. 
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The term 7/ ln(. . . ) in Rp is cancelled by B(Qp) and there is in fact no dependence 
on ftp nor 5f in da^^(Qi, Qj) any more. The IR cancellation is now assured by the term 



2 a - In^-lllns, 



which is implicitly present in B(Qi,qi,q2). 

The situation is, however, not as good as we described above. There is one important 
complication due to the use of the weighted events, at the level of the crude distribution. 
Let us go back again to the case of EEX 




W&c da™ (108 



Now the problem is that photons in 5Q cannot be "hidden" , because WpsR does not obey 
the "IR-safeness" condition 



f 



'11 



lim Wp^(n', k[, . . . , kl, . . . , k n i) — Wp^(n 1, k[, . . . , k\_ x , k' i+1 , . . . , k' n ,) 



Even the softest photons contribute the finite ratio (///), and this contribution is essential 
for the IR-cancellations and for the overall normalization. 

There is, however, a way of saving our method of replacing Qp with Qj. Let us 
repeat again the calculation of eq. (|103|) assuming that photons "hidden" inside 5fl do 
not contribute the factor (///) to the overall weight. We are able to carry out the 
calculation as before, obtaining the modified exponential factor 

00 „ 1 n' 

y{A} = Y, / dr n/+2 (X;k[,...,k' n/ ,q 1 ,q 2 ) — ^25,(^6(0/,^ 

n'=0 H ■ 1=1 (log) 

exp ( L if §f{k) m^/s Q ) ) h{k '^ • • • ' q ^ 

It is very important that the effect due to the omission of (///) in the overall weight can 
be evaluated analytically, and therefore corrected for analytically. In other words we shall 
be able to compensate analytically for the missing average contribution to W^ s r ^ from the 
hidden photons. The evaluation of the integral over Sfl is based on the observation that 



SUk) = S f (k) 



f(0,m}/s) 
mm 2 f /s Q ) 
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where q*, i — 1, 2, are defined such that (g*) 2 = m 2 (sg/s). Furthermore, in the QMS, they 
have the same directions as the original q} and the same total energy, qf + q%° = ^J~sq- 
With the help of the above we get 

r d 3 k ~ f(6,m 2 /s) 
hn = S A k ) f{e J f / SQ) = ^Qf[B(^i, <£, <? 2 *) ~ B{Sl F , <£, g*)]. 

We have at our disposal an analytical representation of the function B(Q, qi, 52), for spher- 
ical Q and regularized with m 7 , in terms of logarithmss and dilogarithms, for arbitrary qi, 
not necessarily antiparallel. Although I$q is IR-finite by definition, it is useful to keep m 1 
and evaluate separately B(Qj) in the CMS and B(Qp) in QMS, where the corresponding 
IR-boundaries are spherical, and subtract the results afterwards. We may calculate B's 
in any frame, because they are (when regularized with m 7 ) essentially Lorentz-invariant 
- we only need to transform IR-boundaries Q correctly from one frame to another. 

Summarizing: in the realistic case of the weighted events (with non-IR-safe weights) 
the method in which we hide photons in 5Q = Qi \ Qf leads to a new crude distribution 
similar to that in eq. ( |106| ) with the new 

Rf = If In ( 5 f SQ+ s 2 Q K ' Q ) + 2Qj« J B(fi / , g*, q* 2 ) - 2Q 2 f aB(tt F , g*, g*). (110) 

As a consequence, the above exponential factor does not cancel out exactly in the model 
weight with the correction to the YFS form factor as before, and we have the following 
additional correcting factor in the model weight: 

Whide = e" 2 "^^ 1 ' 91 *'^^^' 91 *' 92 *^ (111) 

It should be really present in the model weight, but in the program, for historical reasons 
and for convenience, it is included in the crude weight. 

The important profit from the above method is that with the above fix we can now 
make our calculation for the CEEX model with the ISR-FSR interference (IFI) switched 
on. 

Note that the above treatment is more elaborate than the analogous one in BHLUMI 
because it is valid for finite rrif, while in BHLUMI we use the approximation rrif ^/s. 
The correcting weight in BHLUMI is a simple one-line expression while here it expressed 
by a long series of logarithms and dilogarithms. 

Let us finally add a side remark: another valid method of realizing a hypothetical 
Monte Carlo with common Q for ISR and FSR is to generate photons using small Qf 
and to apply a brute-force rejection of all events with one or more photons falling into 
Qj. This was used in the early version of BHLUMI. We do not like this method because 
it may lead to an excessive number of events with zero weight, lowering substantially the 
efficiency of the MC. 
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3.9 Photon multiplicity enrichment 

Let us finally describe yet another complication of the primary distribution, which is 
introduced for technical reasons, i.e. in order to get better total weight distribution, and 
a smaller rejection rate in the process of turning weighted events into unweighted events. 
This modification is not necessary for weighted events. (It can be switched off by adjusting 
input data). 

The problem is essentially due to the introduction of the ISR-FSR interference, which 
we have already called IFI. As already known from the 0(a l ) case, the weight that 
introduces IFI is sharply peaked around 1 and has a strict upper bound Wifj < 2. The 
destructive interference with the weights Wxfx ~ occurs in the backward scattering 
(in the fermion scattering angle) and there is a little constructive interference Wwi ~ 2, 
mostly in the forward direction. The maximum weight = 2 or a factor of 2 in the primary 
primary cross section solves the problem, at the expense of the factor of 2 rejection 



rate 14,31 



In the case of n photons, however, the same leads to: the maximum weight = 2 n or 
equivalently the increase of the primary cross section by a factor of 2 n , and consequently 
increase of the generated photon multiplicity by a factor of 2. Of course, almost all of this 
increase is artificial, and it is compensated by Wifi, leaving only some small net effect 
due to IFI in the MC events. A more sophisticated method would be to increase the 
primary primary cross section (and photon multiplicity) more selectively, that is in a way 
dependent on the fermion scattering angle. For the moment we do not do it. It may be 
done in the future. What we do is the following: we increase the primary cross section 
by a factor of 2 A , where A is not equal to 2, but is adjusted empirically so that the tail of 
the total MC weight is acceptable. We have found that the value A ~ 1.25 is the optimal 
one. 

The introduction of the above A factor affects all formulas for the primary cross section 
and for W in a rather trivial way, so we do not write it explicitly. The only non-trivial 
modification is in the compensating factor for the hidden photons in the previous section. 
This can be understood and implemented as a modification of the electric charge of the 
final fermion Qf — > XQ f in the primary cross section. 

3.10 Entire MC algorithm top-to-bottom 

Specializing to the CEEX model, we can summarize the results of the three previous 
sections as follows 

OO OO „ 

4Lx{A} = £ £ £ / a w^ EX wg™ w hide dogjjcn,). (112) 

f=fi,r,d,u,s,c,b n=0 n'=0 

Here daJ^^iVti) is obtained from the product of the ISR and FSR primary differential 
distributions 

daf^Qj, n F ) = d? p n f \nj) dffi(n F ), (113) 
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see eqs. (|75|) and (j92|) , by means of hiding/ ignoring FSR photons in Consequently, 

(r) 

in the evaluation of Wceex an d of all other weights, only momenta outside the common 
IR-domain Qi enter. 

The value of the integrated cross section with the acceptance function A is obtained 
in the MC run in a standard way 

4Lx{A} = (A W£> EX W hide )a Pn \ (114) 

The acceptance function A may for instance define the entire cross section (A = 1), or 
just a single bin in the histogram of cos 8 for the outgoing fermion, or any other IR-safe 
observable. The overall normalization is based on 

^ = E EE = E / d ^i) E / * ri (^) 

f=^,...,b n=0n'=CT n=0 ^ n'=0 ^ 

1 

CO p p 

E E/^ ri (^)= E /«fo<4 ll («(i-«))J r o(«)7.« vl ^ ) 



f=fj,,...,b n=0 f=H,-,b 







(115) 



where we have exploited the property J ^ d# Pn (fiF) = 1 of eq. (|77|), and the ISR part is 
taken from eq. (OH) . Note that we have put t> max = 1, understanding that crg orn (s) = 



below the threshold, s < Am'j. 

We shall now describe the entire generation of the MC event according to daJ^AQi) 
from the top to the bottom, as is done in the program, starting from the generation of v 
describing the total energy loss due to ISR, the type of final fermion / and the photon 
multiplicities n and n' . Generation of photon energies and angles comes later, using 
methods already described in detail in the previous section. 



3.11 ISR spectrum and fermion type 

First comes the important practical question: Shall we (a) generate first the fermion type 
/ and later v or (b) vice versa? Both options are technically realizable. In case (a) we 
would calculate numerically 

i 

r Pri* _ / J„ _/ ( (i _.,\\ <T( n ,\ => , 1 Te-l ir 7e-7e 



erf 1 = J dv < rn (s(l - v)) J (v) 7 e T/*-V*-* (116) 



for f = fi,r, d, u, s, c, b, and generate fermion type / first; and later on, for a given /, we 
would generate the variable v according to 

da Pri * 

= - «)) Mv) 7et;*-V^. (117) 
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This looks as a natural solution; however, the generation and integration of the ISR 
spectrum daj n /dv is done numerically in a MC module that creates a look-up matrix, 
which memorizes very precisely the shape of the distribution during the initialization 
phase of the MC run (before MC event generation). In this method we would need 
several initializations, creating several tables of this kind. This is feasible, but not very 
convenient. The situation is much worse, when beamstrahlung is switched on because 
in this case the 1-dimensional problem of the generation of v is replaced with the 3- 
dimensional problem of generating v,z\,z 2 , and consequently we would need to manage 
several sets of 3-dimensional look-up matrices. This would make the initialization phase 
rather long in CPU time, and the tables would occupy a lot of processor memory. 

We think that the above scenario is still technically realizable, even in the presence 
of beamstrahlung. Nevertheless, we decided for option (b), which is in our opinion more 
economical. In this case, we generate first the v variable (the case of beamstrahlung is 
described below) according to a distribution summed up over final-state flavour: 

^= £ 4o m (*(l-^)) (118) 

f=H,...,b 

next, for a given v, we generate the final-state flavour / according to a probability 



^BornOU -V)) 

E ^BorrX 1 - V))' 
g=H,...,b 



P f = " Bom P ~ U)) rr . (H9) 



3.12 Inclusion of beamstrahlung 

In the presence of beamstrahlung, the flavour-summed three-dimensional distribution to 
be generated in the very beginning of the MC algorithm is 



Pri* 



dvdz\dz 2 



<h '" 4omW-^^2)Jo(v)%v^- 1 e^V(z u z 2 , v r s y i (i 20 ) 

f=H,...,b 

see also eq. (l26f) . The above 3-dimensional distribution is explored and memorized in the 
look-up matrices in the initialization phase of the MC run. This allows us to generate 
v, Z\, z 2 in a very efficient way for arbitrary 0"B O rn and arbitrary beamstrahlung structure 
function T>(zx, z 2 , \fs). As discussed previously, we admit in V(zi, z 2 , y/s) 5-like singu- 
larities in Zi and, as a result, the MC integration of do Vxi * j dvdz\dz 2 is split into three 
branches with three separate look-up matrices. The above organization assures that the 
beamstrahlung structure function can be a completely arbitrary "user function" . 

Next, for a given set (v, Zi, z 2 ), the final fermion flavour / is generated with the 
probability 

D ._ ^BorX 1 - V)Z!Z 2 ) ^ 12 ^ 



E ^BorrX 1 - V ) Z l Z 2 
=fi,...,b 
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3.13 Photon multiplicities and momenta 

Having denned the total ISR loss variable v or vziz 2 , the ISR photon momenta are gener- 
ated first, and the FSR photon momenta are generated second. For v < e the ISR photon 
multiplicity is zero and for v > e it is generated according to the shifted Poisson distri- 



bution P n -\ with (n) = 7 e ln(u/e), see eq. (|97D , where, in the presence of beamstrahlung 



the modified s' = sz\z 2 enters instead of s into the definition of %, see eq. (0). Next, 
all ISR photons are generated according to the distribution cD^ n of eq. (|92|), with the 
methods already described. The crude weight of eq. ( |93|) is calculated. 

Having generated ISR photons (and optionally beamstrahlung) we now know the total 
four-momentum of the final fermions plus FSR momenta 

n 

X = p a zi + p b z 2 - kj, (122) 

but to start generation of the FSR momenta we need to know only sx = X 2 . First, 
the FSR photon multiplicity n' is generated according to a Poisson distribution with the 
average defined in eq. (|78|). Then, FSR momenta are generated in the rest frame of 
Q = qi + q 2 (QMS). More precisely, their dimensionless energy parameters and angles 
are generated, according to the corresponding FSR primary distribution of eq. (|75|), such 
that the sq = Q 2 needs not be known. The sq is determined as a fraction of sx with the 
help of eq. (pB|), such that photon four-momenta can be constructed in absolute (GeV) 
units in QMS. 

At this point we need to generate angles ip and u in the transformation form QMS 
down to CMS defined in eq. (|59D . Knowing the momentum X = Q — J2j kj in the QMS 
we may apply this transformation and calculate final fermion momenta in the CMS where 
Pi + f>2 = Pa + Pb = 0. In the case of beamstrahlung this transformation bring us to a 
frame where piZi + P2Z2 = 0, and we need an additional boost along beams to brings 
generated momenta to the laboratory system. The same boost is done for ISR photons. 

Removal (hiding) of the FSR photons in 5fl is done at the end of the generation of 
the FSR photons. All remaining photons have the common IR-domain fli defined in the 
CMS where pi +p 2 = p a + Pb = or in the presence of beamstrahlung in the frame where 
PiZi +P2Z2 = 0. 
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4 Structure of the program 



In the following we shall describe the topography of the distribution directory, then the 
programming rules which we follow, and finally we shall briefly describe the functionality 
of the principal modules of the program. 

4.1 Topography of the distribution 

The program source code is organized into modules, also called pseudo-classes, which are 
located in several Unix-type subdirectories of the distribution directory KK-all. The dis- 
tribution directory also contains one additional subdirectory, ffbench, with demonstration 
(template) programs and one subdirectory, dok, with the documentation. The essential 
part of the source code of the /C/C Monte Carlo event generator is located in the two 
subdirectories KK2f and bornv. The tool box of various utilities is located in the subdi- 
rectory glibk, the electroweak library is located in dizet, the MC library of the r lepton 
decays TAUOLA is in the subdirectory tauola, PHOTOS in the subdirectory photos and 
the hadronization package JETSET is in the subdirectory jetset. 



KK2f-all 



KK2f 



bornv 



glibk 



dizet 



tauola 



photos 



jetset 



dok 



ffbench 



Event generator: KK2f, KarLud, KarFin, GPS, 
QED3, HepEvt, PseuMar, TauPair 



Event generator: BornV, MBrA, MBrB, Veskl, 
BStra, FoamA, FoamB, FoamC, VegasA,B,C 



Library of utilities: GLK, MathLib 



Library electroweak corrections: DZface, DIZET 6.21 
sub -directory: tabtest 



Tau-lepton decays: TAUOLA 2.6 



Radiative corrections in decays: PHOTOS 2.02 



Hadronization of quarks: JETSET 7.2 



Documentation in postscript: KKcpc.ps.gz 



Demo programs, example input data, benchmark outputs 
sub-directories: demo, Inclusive, Mu, Tau, Down, Up, 
Botom, Beast 



Figure 1: Topography of the distribution directory 
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4.2 Programming rules 



The program is written in Fortran77 with popular extensions such as long variable names, 
long source lines, etc., which are available on all platforms. In the Makefile in the main 
directory KK-all, there is a collection, of f77 compilation flags, for Linux, AIX, HPUX 
and ALPHA compilers, which should be used to activate these extensions. The program 
is written in such a way that its translation to an object-oriented language such as C++ 
should not be very difficult. In fact the program is divided into modules, which have 
the structure of the C++ classes, as far as it is possible to do it within f77. Below we 
characterize the rules according to which the program was written. 

Each pseudo-class with the name Module consists of a separate source file Module. f 
and the header file Module. h. Each module obeys the following rules: 

• There is only one common block /c_Module/ which contains all class member 
variables, which is placed in the header file Module. h. Each subroutine in the 
Module. f source file includes an INCLUDE ' Module. h' statement. The outside pro- 
grams should never include directly /c_Module/. All input/output communication 
is done with the help of dedicated, easy to use, subroutines. 

• Variables in /c_Module/ are class members and all have the special prefix "m_" in 
their name, for example m_Iterat is the number of iterations. 

• The user has access to some class members through "getters" and "setters"; see 
below. 

• Strong typing is imposed with the help of IMPLICIT NONE. 

• Initializator with the name Module_Initialize performs initialization. Typically 
it initializes variables in /c_Module/. 

• Finalizer with the name Module_Finalize, summarizes the whole run, sets output 
values in /c_Module/, prints output, etc. 

• Maker with the name Module_MakeSomething, or a similar one, does the essential 
part of the job, for instance a maker Module_MakeEvent generates a single MC event. 

• Setter with the name Module_SetVariable is called from the outside world to set 
m_Variable in /c_Module/. For example CALL BornV_SetCMSene(100dO) sets the 
variable m_CMSene=100dO. Only certain privileged variables have a right to be served 
by their own setter, the other ones are in principle "private" . 

• Getter with the name Module_GetVariable is called from the outside world to get 
the m_Variable from /c_Module/. It is a preferred way of sending output informa- 
tion to the outside world. For example, with the help of CALL KK2f _GetXsecMC (xSecPb , 
xErrPb) one gets the MC cross section xSecPb and its error xErrPb in the user pro- 
gram. 

In the following we shall describe all pseudo-classes and their role. 
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4.3 KK2f: Top-level class 



The main purpose of this top-level pseudo-class is to provide the user interface, see section 
^ on the usage of the program. Let us list and explain the main entries in this class: 

• KK2f_ReaDataX('data_file',iReset,imax,xpar) reads the input data file. It should be 
called twice, once with INTEGER iReset=l, for the default data file .KK2f_defaults, 
which is placed in the main directory KK-all. The user has to provide a link to 
this file, or absolute path in the name of the file. (Copying it to the local directory 
is not recommended.) It should be called for the second time for the user data 
./user_data with INTEGER iReset=0, in order to modify some entries in the input 
of the program. The variable INTEGER imax is the dimension of the DOUBLE 
PRECISION xpar(imax). For the moment imax< 3000 is required, but we reserve 
imax= 10000 for future use. The flag INTEGER iReset=0 is for reading data with 
resetting all undefined values to zero, while with iReset=l only entries listed in the 
data file are modified. 

• KK2f_lnitialize(xpar) initializes the whole MC generator. This initializer calls initial- 
izers of other classes like BornVJnitialize, KarLudJnitialize, KarFinJnitialize, QED3_lnitialize, 
GPSJnitialize, TauPairJnitialize. It initializes also the database for the branching over 
final fermion flavours in the class MBrA. Note that BornVJnitialize reads from the 
disk look-up tables for the electroweak corrections and KarLudJnitialize manages 
initialization of the ISR energy spectrum either with the help of Veskl class or, 

in the presence of the beamstrahlung, with the help of the Bstra class, which in 
turn initializes three copies of the Foam package, FoamA, FoamB, FoamC (or of 
Vegas package, VegasA, VegasB, VegasC). 

• KK2LMake generates a single event. It calls KarLudJVIake to make ISR photons, 
KarFin JVIake to make FSR photons and KK2f JVIerge to merge all photons in a single 
list. It invokes KK2f JVIakePhelRand to generate the photon helicities randomly, 
calculates the EEX model weight using CALL QED3JVIake and/or the CEEX model 
weight using CALL GPSJVIake. Optional rejection is performed to produce weight-1 
events and the weight book-keeping is done separately for each final fermion type 
using MBrAJ-ill. Finally, quarks are hadronized using HepEvt JHadronize (interface 
to JETSET), or r decays are simulated with all spin effects (including all spin 
correlations). This is done using subprograms of the TauPair class (interface to 
TAUOLA). 

• KK2f_GetPhotAII(Nphot,PhoAII) provides the user with the momenta of all photons: 
DOUBLE PRECISION PhoAII(100,4) and photon multiplicity INTEGER Nphot. Al- 
ternatively, Nphot is provided by KK2f_GetNphot( Nphot) and the i-th photon mo- 
mentum by KK2f_GetPhotonl(iPhot,Phot), with DOUBLE PRECISION Phot(4). 

• KK2f_GetFermions(ql,q2) provides the user with the momenta of the final fermions 
DOUBLE PRECISION ql(4),q2(4). 
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• KK2f_GetBeams(pl,p2) provides the user with the momenta of the beams DOUBLE 
PRECISION pl(4),p2(4). 

• KK2f_GetWtAII(WtMain,WtCrud,WtSet) can be used to get access to the main MC 
weight WtMain and the list of all alternative weights WtSet(lOOO). The weight for 
the crude differential cross section WtCrud is also provided. All of them are DOUBLE 
PRECISION type. Alternatively, the getter KK2f_GetWt(WtMain,WtCrud) may be 

more convenient. 

• KK2f_Finalize may be called at the end of the MC run, in order to perform the final 
book-keeping and printing. 

• KK2f_GetXSecMC(XSecPb,XErrPb) should be called after calling KK2f_Finalize in 

order to get the total cross section (in picobarns) and its absolute error: DOUBLE 
PRECISION XSecPb, XErrPb. 

• KK2f_GetVersion(Version) and KK2f_GetDate(Date) provide the user with the version 
number DOUBLE PRECISION Version and the release data CHARACTER*14 Date 

information. This should help to keep track of the development of the program. 

• KK2f_Print(iel,ie2) can be used to print the actual MC events, limiting their serial 
number to stay between INTEGER iel and INTEGER ie2. 

There are several other getters in the KK2f class, which are mainly for internal use. 

Let us briefly list other subroutines in the KK2f class, which are not called by the 
user of the program: 

• KK2f_Wignerlni(KFbeam,CMSene,PolBeaml,PolBeam2, Polarl,Polar2) does Wigner ro- 
tation for spin polarization vectors of beams. Beam polarization vectors (in input 
data) are defined in the beam particle rest frames, which are reached from the CMS 
by a simple z-boost without any rotation. (The first beam is parallel to the z-axis.) 

• KK2f_ZBoostAII(exe) performs a z-boost on all momenta of the event. This z-boost 
corresponds to beamstrahlung or beam spread and is done at the very end of gen- 
eration, after the calculation of the matrix element. 

• KK2f_Dsig0verDtau(mout,Rho) is only for documentation and testing purposes. It 
calculates the distribution da/dr corresponding to WtCrud, normalized with respect 
to dr = Lorentz invariant phase space. 

• KK2f_Merge merges lists of ISR and FSR photon momenta. The resulting merged 
photons are ordered according to their energy. 

• KK2f_MakePhelRand generates the photon helicities randomly. 
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4.4 HepEvt: HEP event record class 



The pseudo-class HepEvt has the double purpose of (a) being another user interface, 
alternative to getters in KK2f, and (b) managing also the hadronization of quarks. The 
user may also traditionally put the common block of the HepEvt class directly into its 
code. It is a DOUBLE PRECISION version of the standard PDG /HEPEVT/ common 
block, for a maximum of 2000 particles. 



Let us now list and explain the subprograms in the HepEvt class. 

• HepEvt_Fill fills in all of common block /cLHepEvt/ using the subprogram Hep- 
Evt_Fill. Momenta are provided by getters from KarLud and KarFin. 

• HepEvt_Fill writes a single particle record into the common block /cLHepEvt/. 

• HepEvt_Hadronize(HadMin) arranges jets with the help of LuJoin and LuShow, and 

hadronizes quarks using LuExec. 



• CALL HepEvt_GetBeams(pl,p2) provides four-momenta of the two beams DOUBLE 
PRECISION pl(4),p2(4). 

• CALL HepEvt_GetFfins(ql,q2) provides four-momenta of the two final fermions DOU- 
BLE PRECISION ql(4),q2(4). 

• CALL HepEvt_GetPhotAII(NphAII,PhoAII) provides the photon multiplicity INTEGER 
NphAII and the photon four-momenta DOUBLE PRECISION PhoAII(100,4). 

• HepEvt_GetPhotBst(nPhot,Phot) provides the multiplicity INTEGER nPhot (=0,1,2) 
and four-momenta DOUBLE PRECISION Phot(100,4) of the beamstrahlung photons. 

• HepEvt_LuHepc This is the double-precision version of LUHEPC of JETSET. It trans- 
lates DOUBLE PRECISION /c_HepEvt/ into the old style REAL*4 Lund commons. 



* 



INTEGER nmxhep ! maximum number of particles 

PARAMETER (nmxhep=2000) 
DOUBLE PRECISION phep, vhep 

INTEGER nevhep, nhep, isthep, idhep, jmohep, jdahep 
COMMON /d_HepEvt/ 



$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 



vhep (4, nmxhep) ! vertex [mm] 
SAVE /d_hepevt/ 



nevhep, ! serial number 

nhep, ! number of particles 

isthep (nmxhep) , ! status code 

idhep (nmxhep) , ! particle ident KF 

jmohep (2 , nmxhep) , ! parent particles 

j dahep (2 , nmxhep) , ! childreen particles 

phep (5 , nmxhep) , ! four-momentum, mass [GeV] 



* 
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4.5 KarLud: Crude level MC for ISR 



The simulation of ISR together with the beamstrahlung and generation of the type of 
final-fermion flavour is implemented in the KarLud pseudo-class. 
Let us now list and explain all subprograms in the KarLud class: 

• KarLud_lnitialize(xpar_input,XCrude) is the initializer of the class. It initializes the 
generation of v with the help of VesklJnitialize or of v,zi,Z2 with the help of 
BStraJ initialize. It defines the primary integrated cross section on which the en- 
tire normalization is based. 

• KarLucLSmearBeams implements beam spread. This is correct only for a small spread 
< 2 GeV. It should not be used together with beamstrahlung, since this has not been 
tested yet. The distribution is Gaussian p(X) = N exp((X-CMSene/2) 2 /(2DelEne 2 2)) 
(DelEne is the dispersion of the beam energy Ebeam, not of CMSene). 

• Karl_ud_Make(PX,wt_ISR) generates ISR photons with the help of other subprograms; 
PX(4) is the four-momentum left after photon emission, wtJSR is the ISR component 
of the "crude weight". The other subroutines called here are: Karl_ud_SmearBeams 
(see below), BornV_SetCMSene(XXXene), which resets the CMS four-momentum in 
BornV in the case of beam smearing, Veskl_Make, which generates the variable v, 
in the absence of beamstrahlung (alternatively it is done with VegasA_Generate for 
KeyFix=2 ), or BStra_Make, which generates v,z\,Z2 in case of beamstrahlung. The 
ISR photons are generated with the help of Karl_ud_YFSini, see below, and the type 
of final-state fermion KF is generated with the help of MBrA_GenKF. Finally, if FSR 
is off, then final momenta are produced locally with the help of Kinl_ib_phspc2 

• Karl_ud_Finalize(Mode, XKarlud, KError) calculates the crude cross section XKarlud 
and its error KError, and prints out final statistics. The crude cross section is 
coming from the Veskl_Finalize (alternatively from VegasA_GetlntCrude) or, in case 
of beamstrahlung, from BStra_GetXCrude, 

• Karl_ud_YFSini(XXXene,vv, PX.Wtlni) is generating the ISR photon momenta. Its 
input is the total energy available XXXene, and vv= v. Wtlni is the ISR com- 
ponent of the "crude weight" and PX(4) is the four-momentum left after photon 
emission. Karl_ud_YFSini calls BornV_GetAvMult to get the average ISR multiplicity, 
Karl_ud_PoissGen and Karl_ud_AngBre, see below. 

• Karl_ud_PoissGen generates photon multiplicity. 

• Karl_ud_AngBre generates photon angle. 

• Karl_ud_ZBoostAII performs z-boosts of all photons. 

• KarLud_GetPhotons(nphot,sphot) provides all ISR photons. 

• Karl_ud_GetPhotonl(iphot,phot) provides single ISR photons. 
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• KarLucLGetPX(PX) provides four-momentum PX, see above. 

• Karl_ud_GetBeams(pl,p2) provides the beam momenta. In the case of beamstrahlung 
pl,p2) are beams after beamstrahlung. 

• Karl_ud_GetBeasts(pl,p2) provides the collinear photons of the beamstrahlung. 
4.6 KarFin: Crude level MC for FSR 

The FSR is implemented in the KarFIN class. This package was already used in the 
KORALZ |]J program for some time. Thanks to recent improvements it now works 
properly without any approximations close to r threshold. The rrif ^/i approximation 
is not used anymore. 

Let us now list and explain all subprograms in the KarFin class: 

• KarFinJnitialize initializes some internal variables (weight book-keeping). 

• KarFin_Make(PX,amfil,amfi2,CharSq,WtFin) generates FSR photons with the help 
of other subprograms, see below. PX is the four-momentum of the entire FSR 
system (fermions + photons), amfil,amfi2 are masses of the final charged pair (not 
necessarily equal), CharSq is the final-state fermion charge squared and WtFin is the 
FSR part of the crude weight. 

• KarFin_YFSfin(PX,amfil,amfi2, CharSq, WtFin) generates momenta of the FSR pho- 
tons. PX, amfil, amfi2, CharSq, WtFin are as defined above. It calls KarFin_PoissGen 
and KarFin_AngBre to generate photon multiplicity and angles, then KarFin_Kinfl 
and and KarFin_Piatek see below. 

• KarFin_Kinfl(PX,...,phsu) transforms from the rest frame of Q=ql+q2 QMS down 
to the laboratory through the intermediate rest frame of PX=ql+q2+ phsu. 

• KarFin_Piatek( Masl,Mas2,CharSq,WtMlist, Wt3)[] optionally removes photons be- 
low i? m in from the list of photons, appropriately modifying the crude weight. Masl,2 
= fermion masses, WtMlist = list of mass weights (///) for all photons. Wt3 = prod- 
uct of (///) for the alive (not hidden) photons. The correcting weight is calculated 
with the help of BVR_Btildc, calculating B for q* and BVR_Btilda for q { . 

• KarFin_PoissGen generates the photon multiplicity randomly. 

• KarFin_AngBre generates the photon angles randomly. 

• KarFin_Kinfl transforms the final fermions and all photons from QMS through the 
Z-frame to CMS. Random rotation with angles is applied in the intermediate 
rest frame of PX [Z boson) using KarFin_BostEul. 

21 Written in CERN, Piatek=Friday, 22.IX.1989 (S.J.) 
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• KarFin_BostEul performs Lorentz transformations consisting of: (1) parallel boost 
from the final fermions rest frame to fermions + photons rest frame (Z frame); (2) 
two rotations with angles ip,u;; (3) parallel boost to the laboratory system CMS. 

• KarFin_ZBoostAII(exe) does an additional z-boost of all particles in case of beam- 
strahlung. 

• KarFin_Finalize prints final statistics. 

• KarFin_GetNphot(nphot) provides the FSR photon multiplicity. 

• KarFin_GetPhotonl(iphot,phot) provides the four-momentum of a single FSR photon. 

• KarFin_GetPhotons(nphot,sphot) provides the four-momenta of all FSR photons. 

• KarFin_GetFermions(qfl,qf2) provides the four-momenta of the final fermions. 

• KarFin_WtMass(WtMass) provides the product of (///) for the alive (not hidden) 



4.7 BornV class: particle data base and ISR spectrum 

Class BornV is serving as a data-base for fermion properties such as mass, charge, isospin, 
colour and other fermion-type dependent parameters relevant to MC generation, like the 
maximum weight for rejection. It also reads from the disk and keeps the EW form factors 
produced by the interface to DIZET 6.21. 

All other classes use the data-base of the BornV class through its getters, see be- 
low. The data-base is located in the class common block c_BornV, which is initialized by 
BornVJnitialize from the default input data file KK-all/.KK2f_defaults passed by arguments 
from KK2f_lnitialize. Optionally, only after user modifications, see section ^ on the usage 
of the program. 

For this particular class it is instructive to look into the list of the class member 
variables in the class common block/c_BornV/. Below we quote part of the BornV. h 
source code: 



photons. 



* 



COMMON /c_BornV/ 



* Tables of EW formf actors 

$ m_cyy (m_poinl+l ,7 , 16) , 

$ m_czz(m_poin2+l ,7 , 16) , 

$ m_ctt(m_poin3+l,m_poinT+l,7,16) , 

$ m_clc(m_poin4+l,m_poinT+l,7,16) , 

$ m_syy (m_poinl+l , 16) , 

$ m_szz(m_poin2+l , 16) , 

$ m_stt(m_poin3+l,m_poinT+l,16) , 

$ m_slc(m_poin3+l ,m_poinT+l , 16) , 

$ m_GSW(100), ! form-factors, 



at the actual energy/angle 



formf actor, table 
formf actor, table 
formf actor, table 
formf actor, table 
QCD correction, table 
QCD correction, table 
QCD correction, table 
QCD correction, table 
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$ m_QCDcor, 



! QCD correction, at the actual energy/angle 
EVENT 





in 


CMSene 


1 Tnitial valnp of ("IMS c*T)C*ttct\t 

. ±111 L> _L CI J- V CJ.. _L L1C Ul 1U did tl y 


$ 


m. 


.XXXene , 


! CMS energy after beamstrahlung or beam spread 


$ 


m_ 


xl, 


! 1-zl = xl for first beam(strahlung) 


$ 


m_ 


_x2, 


! l-z2 = x2 for second beam(strahlung) 


$ 


m_ 


_vv, 


! v = 1-sprim/s 


$ 


m_ 


.AvMult , 


! Average photon multiplicity CRude at given v 


$ 


m_ 


.YFSkon, 


! YFS formfactor finite part 


$ 


m_ 


_YFS_IR, 


! YFS formfactor IR part 


$ 


m_ 


.vvmin, 


! minimum v, infrared cut 


$ 


m_ 


.vvmax , 


! maximum v 



$ m_HadMin, 

* Basic QED 

$ m_alfinv, 

$ m_alfpi, 

$ m_Xenph , 

* EW parameters 

$ m_MZ , 

$ m_amh, 

$ m_amtop, 

$ m_swsq, 

$ m_gammz , 

$ m_amw, 

$ m_gammw , 

$ m_Gmu , 



minimum hadronization mass [GeV] 



1/alphaQED, Thomson limit 
alphaQED/pi 

Enhancement factor for Crude photon multipl. 

Z mass 
Higgs mass 
Top mass 
sin(thetaW)**2 
Z width 
W mass 
W width 

Fermi constant (from muon decay) 



* Table of fermion parameters, quarks (l->6) and leptons (11->16) 



$ 
$ 
$ 
$ 
$ 
$ 
$ 



m_KFf erm(20) , 
m_NCf (20) , 
m_Qf (20) , 
m_T3f (20) , 
m_helic(20) , 
m_amf erm(20) , 



fermion KFcode (l->6) and (11->16) 
number of colours 
electric charge 
isospin, L-hand component 
helicity or polarization 
fermion mass 
auxiliary parameter 

! Generation flag, only for SemiAn. ! 



m_auxpar(20) , 
$ m_IsGenerated(20) , 

* Normalisation 

$ m_gnanob, ! GeV~(-2) to nanobarns 

* Initial/final fermion types 



$ m_KFini , 

* Test switches 

$ m_KeyINT, 

$ m_KeyElw, 

$ m_KeyZet , 

$ m_KeyWtm, 

$ m_out 



! KF code of beam 

ISR/FSR interef erence switch 
Type of Electroweak Library 
Z-boson on/off 

Photon emission without mass terms 
output unit for printouts $ 
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us now list and explain all subprograms in the BornV class: 
» BornVJnitialize(xpar) initializes data members in /c_BornV/. 

» BornV_StartEW(xpar) initializes electroweak formfactors in /c_BornV/. There are 
two versions of this routine: one in KK-all/bornv/BornV_StartEW.f which reads elec- 
troweak formfactors from the disk file and another one in KK-all/dizet/BornV_StartEW.f 
which calculates it using DIZET library. See section 4.13 for more details. 

) Born V_Read All reads from the disk-file pretabulated EW form factors for /i and r 
leptons, and for d, u, b quarks. For s, c quarks the form factors of d, u are used. 

» BornV_ReadFile(DiskFile,KFfin) reads from the disk a single file for a single final 
fermion. 

» BornV_StartDZ(xpar) Initialized DIZET library using current input data in xpar. 

» BornV_ReBinl, BornV_ReBinla, Born V_ReBin2, BornV_ReBin2a subroutines map the 
variable r G (0,1) (random number) into v G (0,t> max ). Various methods are used 
to do it, with various kinds of the mapping function. 

. DOUBLE PRECISION FUNCTION BornV_RhoFoamC(xarg) is the integrand for FoamC 
in the 3-dimensional mode for beamstrahlung. Remember that BornV_Crude and 
BornV_MakeRho use the hidden input m_XXXene. BornV_Crude is in the R-units 
(point-like cross-section at v /i=m_XXXene). It defines m_vv, which is later on ex- 
ported to KarLud. 

• DOUBLE PRECISION FUNCTION BornV_RhoFoamB(xarg) is the integrand for FoamB 
in the two-dimensional mode for beamstrahlung (it defines m_vv). 

. DOUBLE PRECISION FUNCTION BornV_RhoFoamA(xarg) is the integrand for FoamA 
in the one-dimensional mode for beamstrahlung off and on (it defines m_vv). 

• DOUBLE PRECISION FUNCTION BornV_RhoVeskol(R) is the integrand of Veskol. 
(The comment about hidden input m_XXXene applies.) In the case of beamstrahlung 
the additional normalization factor Circee(ld0,ld0) is added inside BStraJnitialize 
(it defines m_vv). 

» BornV_MakeGami(CMSene,gamiCR, garni) calculates GamiCR= 7 e and gami= 7 e as 
functions of CMSene. 

» Born V_MakelSR(Rho) This procedure is tightly related to ISR photon generation in 
KarLud. It provides Rho(m_vv, m_XXXene), the primary distribution of v. It also 
calculates m_AvMult, which is later used in KarLud_YFSini; m_YFSkon ,m_YFS_IR, 
which are later used in GPS_Make and QED3_Make. 
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• DOUBLE PRECISION FUNCTION BornV_Crude(vv) calculates the crude total Born 
cross section summed over fermion types. It exploits the fact that the Born differ- 
ential distribution reads a + bcosO + dcos 2 8. (Hidden input is m_XXXene). It is 
used in Born V_RhoVeskol, BornV_RhoFoamA, etc. 

• DOUBLE PRECISION FUNCTION BornV_Differential(Mode,KFf,svar,CosThe,...) is 

the Born differential distribution. For Mode=0 it is a crude version of pure Born, 
no spin, no EW corrections. For Mode=l it is the full result with EW corrections 
spin, etc. In this mode it is used in QED3, and for all kinds of tests. For Mode=3 
it is used in the tests of pretabulation. In this case GSW(s,theta) has to be provided 
from the outside, with the help of BornV_SetGSW0. 

• DOUBLE PRECISION FUNCTION BornV_Simple(KFi,KFf,svar,costhe) provides for 
BornV_Differential an unsophisticated Born differential distribution without EW cor- 
rections, with the Z and 7 s-channel exchange. 

• DOUBLE PRECISION FUNCTION BornV_lntegrated(KFfin,svar) is used only in semi- 
analytical programs. It calculates the total Born cross section. For KFfin = it sums 
over all allowed flavours; otherwise, for KFfin. NE.O, it calculates the cross section 
for the actual value of m_KFfin. 

• DOUBLE PRECISION FUNCTION BornV_Dizet(Mode,KFi,KFf,svar,CosThe,...) pro- 
vides for BornV_Differential the differential Born cross section with/without EW 
corrections. For Mode=0 it provides pure Born and for Mode=l electroweak correc- 
tions are added. KFi.KFf can also be negative for an antiparticle; in this case it is 
important to produce tables with the correct input KFini, KFfin. 

• BornV_lnterpoGSW(KFf,svar,CosThe) calculates EW form factors from look-up ta- 
bles, using linear interpolation. 

• BornV_givizo(idferm, ihelic, sizo3, charge, kolor) provides electric charge, weak isospin 
and colour of the fermion, where idferm =1,2,3,4 denotes: neutrino, lepton, up, down 
quark; negative idferm=-l,-2,-3,-4, denotes the corresponding antiparticle; ihelic = 
+1, —1 denotes right- and left-handedness (chirality) sizo3 is the third projection 
of weak isospin (±1/2), charge is the electric charge (in units of magnitude of the 
electron charge), kolor is the QCD colour, 1 for lepton, 3 for quarks. 

• DOUBLE PRECISION FUNCTION BornV_SigOnb(CMSene) provides the point-like 
muon cross section in nanobarns for the normalization purpose. 

Communication subprograms (setters and getters) used by all other classes are the fol- 
lowing: 

22 Note that in the test mode KeyEwl=0 and Mode=l we use BornV_Simple, which will perhaps have 
to be changed in the future, because of the lack of spin effects. At this stage, however, we are bound to 
use it because the KeyZet etc. are implemented only in BornV_Simple and not in BornV_Dizet. 
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• BornV_GetParticle(KFferm, mass, Qf, T3f, NCf); for the fermion type INTEGER KF- 
ferm provides its QCD colour INTEGER NCf, mass, electric charge and weak isospin 
DOUBLE PRECISION mass, Qf, T3f. 

• DOUBLE PRECISION FUNCTION BornV_GetMass(KFferm); for the fermion type IN- 
TEGER KFferm provides its mass. 

• DOUBLE PRECISION FUNCTION BornV_GetCharge(KFferm); for the fermion type 
INTEGER KFferm provides its electric charge. 

• INTEGER FUNCTION BornV_GetColor(KFferm); for the fermion type INTEGER KF- 
ferm provides its QCD colour. 

• DOUBLE PRECISION FUNCTION BornV_GetAuxPar(KFferm); for the fermion type 
INTEGER KFferm it provides its auxiliary parameter. 

• BornV_SetKeyElw(KeyElw) sets the EW switch KeyElw. 

• BornV_GetKeyElw(KeyElw) gets the EW switch KeyElw. 

• BornV_GetKeyZet(KeyZet) sets the Z boson switch KeyZet. 

• BornV_SetKeyZet( KeyZet) gets the Z boson switch KeyZet. 

• BornV_SetCMSene(CMSene) sets the CMS total energy CMSene. 

• BornV_SetMZ(MZ) sets the Z boson mass. 

• BornV_GetMZ(MZ) gets the Z boson mass. 

• BornV_GetGammZ(GammZ) gets the Z boson width. 

• BornV_GetGmu(Gmu) gets the GFcrmi- 

• BornV_GetSwsq(Swsq) gets the electroweak mixing angle. 

• BornV_GetAlflnv(Alflnv) gets the c^qed/ 71 "- 

• BornV_GetAvMult(AvMult) provides the average ISR multiplicity. 

• BornV_GetYFSkon(YFSkon) provides the finite part of the YFS form factor. Used in 
QED3. 

• BornV_GetYFS_IR(YFS_IR) provides IR (cut-off-dependent) part of the ISR YFS form 
factor. Used in QED3. 

• BornV_GetQCDcor(QCDcor) provides the QCD correction factor, defined by DIZET. 

• BornV_GetVV(vv) provides v=w. 
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• BornV_GetVXX(vv,xl,x2) provides v=vv, X{ = 1 — z\ = xl,x2. 

• BornV_GetGSW(GSW) provides the EW form factors GSW(k) k=l 7. It is used in 

GPS and in BornV_Dizet. Note that BornVJnterpoGSW has to be called before, in 
order to interpolate properly. 

• BornV_SetGSW(GSW). For special tests of pretabulation the values of EW form 
factors GSW(k) k=l,...,7 can be set with this subprogram from outside. 



4.8 Bstra, IRC and MBrB classes for beamstrahlung 

As was already described, the MC integral for beamstrahlung and ISR has three com- 
ponents: J dvdz\dzii J dvdzi and j dv (J dvdz 2 is obtained by symmetrization). The 
corresponding three-fold branching method is managed by the class MBrB. The class 
Bstra contains mainly the interface to the Foam and Vegas packages, which generate in 
each branch the corresponding subset of the variables v, z\,z%. 

Since Foam makes look-up tables for the distribution to be generated, we need three 
copies of the Foam, which are FoamA, FoamB, FoamC. Similarly for Vegas. There are 
only five subprograms in the Bstra class: 

• BStra_lnitialize(KeyGrid,Xcrude) initializes FoamA, FoamB, FoamC (or Vegas A, 
VegasB, VegasC) and the MBrB for book-keeping in the three-fold branching. 
It also provides the primary integrated cross section used to establish the overall 
normalization in Karlud and KK2f. 

• BStra_Make(vv, xl, x2, MCwt) randomly chooses one of the branches with the help 
of MBrB_GenKF. 

• BStra_Finalize(lnteg,Errel) is calculating the total cross section using the average 
provided by MBrB_MgetAve. This is for control only. 

• BStra_GetXCrude(XCrude) provides the primary integrated cross section that enters 
the overall normalization. 

• BStra_GetlntegMC(lntegMC,ErRelMC) provides the value of the primary integrated 
cross section from the entire MC run. This is for control only. 

The CIRCE library of the beamsstrahlung structure functions [f23J is placed in the IRC 
module. The only changes are: adding prefix IRC_ to names of all subprograms and the 
common block /circom/ is renamed as /c_IRC/. It is exploited in the BStra class. 



4.9 BVR class: virtual corrections 

The class BVR is a collection of the complex functions used in the calculation of the vir- 
tual corrections. It provides also the B function necessary to calculate the compensating 



weight in the procedure of the removal of the FSR photons, see subsection pT8| . It has its 
own library of complex logarithms and dilogarithms. 



63 



4.10 QED3 class: EEX distributions 



The EEX differential distributions is implemented in the QED3 module. This module is 
rather monolithic. It contains the initializer QED3_lnitialize, the maker QED3_Make and 
several small functions for calculating virtual corrections and up to third-order leading- 
logarighmic structure functions. The basic ingredient in the EEX differential distributions 
is the Born differential distribution that comes from the BornV class. 

4.11 GPS class: CEEX Matrix element 

The CEEX matrix element is programmed in the GPS class. It calculates spin amplitudes 
for the e~e + — > ffwy process. It has grown to a very large module (almost 5000 lines of 
code) and will therefore be split in the next version into a low-level library of GPS tools 
and the module CEEX, which calculates solely the CEEX spin amplitudes. 

The main subprogram in the GPS class is GPS_Make, which calculates 0(a r ) r = 
0, 1, 2 CEEX spin amplitudes m_AmpExpo0(4,4,4,4), m_AmpExpol(4,4,4,4), m_AmpExpo2(4,4,4,4) 
of the DOUBLE COMPLEX type. Photon helicities are generated randomly in the upper 
class KK2f and provided with the getter KK2f_GetPhel. Virtual corrections are provided 
by subprograms from the BVIR class. Masses, charges, isospin of the particles are pro- 
vided by the getters of the BornV class and the electroweak form factor come from 
DIZET 6.21 through the interface subprogram BornV_GetGSW. The three model weights 
0(a r ), r = 0, 1,2, are calculated in GPS_Make using GPS_MakeRho for polarized beams 
and unpolarized final fermions. Subprogram GPS_MakeRho2 is calculating model weights 
for polarized beams and polarized final fermions. It is used by TaupairJmprintSpin to 
implement spin effects in r decays. 

Let us now list and explain three groups of subprograms in the GPS class, (a) the 
main subprograms calculating the spin amplitudes, (b) the library of basic tools and (c) 
the communication subprograms (setters and getters). The first group includes: 

• GPSJnitialize initializes of the class. It sets some coupling constants, Pauli matrices, 
the axial gauge vector (3. 

• GPS_Make is the main routine that calculates spin amplitudes m_AmpExpo0, m_AmpExpol 
and m_AmpExpo2. Spin amplitudes are calculated in such a way that they are first 

set to zero and then, in the sum over partitions, they are incremented by (3q with 
the help of GPS_BornPlus, by fi[ r) with the help of GPS_HiniPlus and GPS_HfinPlus, 
and by tfp with help of GPSJHiiPlus, GPSJHffPlus and GPSJHifPlus. Three model 
weights are calculated using GPS_MakeRho and set to m_WtSet(i), i=l ,2,3 (or to 
m_WtSet(i),i=51, 52, 53 if ISR-FSR interference is switched off). The best weight 
is set as m_WtBest = m_WtSet(3) (or m_WtSet(53)). Weights are available through 
getter GPS_GetWtSet, see below. 

• GPS_MakeRho(ExpoNorm) calculates differential distributions (normalized to Lorentz- 
invariant phase space) from spin amplitudes m_AmpExpo2, i — 1,2, 3, for polarized 
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beams and unpolarized final fermions. Beam polarizations are set from outside with 
the help of GPS_SetPol Beams. 

• GPS_MakeRho2(wtO,wtl,wt2) is used in TaupairJmprintSpin and it calculates the dif- 
ferential distributions (normalized to Lorentz-invariant phase space) from spin am- 
plitudes m_AmpExpoi, % — 0, 1,2, for polarized beams and 'polarized final fermions. 
Final-state polarimeter vectors are set from outside with the help of GPS_SetHvectors. 

• GPS_BornPlus calculates the spin amplitudes of [3q \ It is optimized for summation 
over partitions. Virtual corrections (boxes and vertices) are included. 

• GPS_Born provides the Born spin amplitudes used in the construction of the hard 
non-IR parts: in GPSJHiniPlus, GPSJHfinPlus and other subprograms. It is essen- 
tially a simplified clone of GPS_BornPlus. 

• GPS_EWFFact creates form factors for electro- weak corrections. They are in vector 
couplings (multiplied by the correcting factors). Because of the cos 6* dependence of 
WW boxes, we had to introduce CosThetD parameter. 

• GPS_HiniPlus calculates the IR-finite part of 1-photon amplitudes for ISR Pi^y- (It 
is equivalent to a testing subprogram GPSJHini). 

• GPS_HfinPlus calculates the IR-finite part of 1-photon amplitudes for FSR P^ y- (It 
is equivalent to a testing subprogram GPSJHfin). 

• GPS_HffPlus calculates the IR-finite part of 2-photon amplitudes for FSR Pi^ 00 y 

• GPS_HiiPlus calculates the IR-finite part of 2-photon amplitudes for ISR /3^ U }. 

• GPS_HifPlus calculates the IR-finite part of 2-photon amplitudes /3^ 10 y for one ISR 
and one FSR photon. 

Let us now list and explain the subprograms that play a role of the library of basic tools. 
(In the future version to be isolated as a separate class, or even several classes). This 
group includes: 

• GPS_PartitionStart(nphot,last) initializes the first partition in the sum over partitions. 

• GPS_PartitionPlus updates the partition vector mJsr, checks if it is the last partition. 

• GPS_BornZero(AmpBorn) sets AmpBorn to zero. 

• GPS_BornCopy(AmpBorn,AmpBorn2) copies AmpBorn into AmpBorn2. 

• GPS_BornSumSq(AmpBorn,Sum) sums up AmpBorn amplitudes squared. 

• GPS_TralorPrepare prepares transformation for Tralor, according to GPS rules. The 
resulting Lorentz transformation matrix is stored for multiple use. 
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GPS_GPS(xi,eta,Rot) defines the basis vectors ei,e2,e3 from £ and r] according to 
GPS rules. Columns in the matrix Rot are ei,e2,es. This subprogram is called in 
GPS_TralorPrepare. 

GPS_TralorDolt(id,pp,q) transforms the four-vector pp from rest frame of fermion id 
to LAB, q is the result. It uses a Lorentz transformation, prepared and memorized 
in the subprogram GPS_TralorPrepare, which has to be called first. This organization 
saves CPU time in the case of multiple calls for several r decay products. 

GPS_TralorUnDo(id,pp,q) is the inverse of GPS_TralorDolt. It transforms pp from the 
laboratory to the rest frame of the final fermion, q is the result. It uses a Lorentz 
transformation, prepared and memorized in the subprogram GPS_TralorPrepare, 
which has to be called first. 

GPS_TraJacobWick(Mode,QQ,pp,rr) is for tests only. It is a Tralor-type transfor- 
mation for the classical Jacob- Wick quantization axes. Not optimized. 

GPS_RmatMake is for tests only. It translates Born spin amplitudes into a double- 
spin density matrix m_AmpBorn — > R ab 

GPS_MakeU(ph, sigma.pl, ml, p2,m2,U) builds the transition matrix U, (u fu). 

GPS_MakeV(ph, sigma.pl, ml, p2,m2,V) builds the transition matrix V, (v fv). 

GPS_MakeUb(ph, sigma.pl, ml, p2,m2,U) builds the transition matrix U, (u fu). 

GPS_MakeVb(ph, sigma.pl, ml, p2,m2,V) builds the transition matrix V, (v fv). 

GPS_MatrU(Cfact,ph, sigma.pl, ml, p2,m2,U) builds the transition matrix U, (u fu). 

GPS_MatrV(Cfact,ph, sigma.pl, ml, p2,m2,V) builds the transition matrix V, (v fv). 

GPS_MatrUb(Cfact, ph.sigma.pl, ml, p2,m2,U) builds the transition matrix U, (u fu). 

GPS_MatrVb(Cfact,ph, sigma.pl, ml, p2,m2,V) builds the transition matrix V, (v fv). 

DOUBLE COMPLEX FUNCTION GPS_Sofl(sigma,ph,pf) calculates the single soft 
photon contribution to the s-factor. 

DOUBLE COMPLEX FUNCTION GPS_Soflb(sigma,ph,pf,mf) calculates the single 
soft photon contribution to the s-factor. 

DOUBLE COMPLEX FUNCTION GPS_soft(sigma,ph,pl,p2) calculates the two-fermion 
s-factor. 

DOUBLE COMPLEX FUNCTION GPS_bfact(sigma,phot,pferm) calculates the diago- 
nal element of the [/-matrix for the massive fermion (the numerator in the s-factor). 
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• DOUBLE COMPLEX FUNCTION GPS_softb(sigma,ph,pl,ml,p2,m2) calculates the s- 
factor. 

• DOUBLE COMPLEX FUNCTION GPS_bfacb(sigma, phot, pferm, mass) calculates the 
diagonal element of the [/-matrix for massive fermion (the numerator in the s- 
factor). 

• DOUBLE COMPLEX FUNCTION GPS_iProdl(L,p,q) calculates the basic inner prod- 
uct of spinors S\(p, q) = U\(p)u^\(q). We exploit the identity q) = — [s + (p, q)}*. 

• DOUBLE COMPLEX FUNCTION GPS_iProd2(Lamp,p,mp,Lamq,q,mq) calculates the 
general spinor product S\ 1 ,x 2 (p,q) for massive spinors u and/or v; mp and mq are 
the masses of four- vectors p and q. Negative mass means an antiparticle. 

• DOUBLE PRECISION FUNCTION GPS_XiProd(p,q) is the auxiliary function called 
in GPS_iProd2. 

The last group includes communication subprograms (setters and getters) and some mis- 
cellaneous routines for debugging: 

• GPS_BPrint(nout,word,AmpBorn) prints 16 spin amplitudes of AmpBorn in a nice 
format on output unit nout. 

• GPS_GetXi(xi,eta) provides £, the basic light-like vector in the laboratory frame, 
entering the definition of all spinors (called k in Kleiss-Stirling papers). 

• GPS_SetKeyArb(KeyArb) GPS_GetKeyArb(KeyArb) sets KeyArb, which is switching 
on/off the use of m_b, KeyArb=0 means f3 — > £. 

• GPS_Setbl GPS_Setb2 switches the axial gauge vector /3=b to another predefined 
value. This is for testing the gauge invariance of the spin amplitudes. 

• GPS_GetWtSet(WtBest,WtSet) provides a complete list of weights. 

• GPS_SetKeylNT(KeylNT) sets the IFI switch KeyINT 

• GPS_SetPolBeams(PolBeaml,PolBeam2) sets the beam polarization vectors. One 
should not forget the Wigner rotation to the GPS frame! 

• GPS_SetHvectors(HvecFerl,HvecFer2) sets the final-fermion polarimeter vectors. 
4.12 TAUOLA and PHOTOS 

TAUOLA and PHOTOS are placed in KK-all/tauola and KK-all/photos. They commu- 
nicate with the rest of the program through an interface class Taupair located in KK- 
all/KK2f. The initialization is performed in the Tauface class as well. The other, very 
important role of Taupair is to implement spin effects in the decays of both r's, including 
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all spin correlations with the rejection method according to the special spin weght; as in 



KORALB [14]. The spin weight is: 



_ a- r ,X A ,X A i,j,l,m 
"spin 



2^ £ 1 £ 2 \ aXa ° \ b \ b JJLn Uoioa • • • «rj JJln \.\<n ' • • «rj 



l_ m U j,m 
a A c A c <7 A d \, /l 3 ft 4 



2 



(123) 



Note that in the present version of the program we include at this step not only the final- 
state spin effects but also beam-polarization effects. This is not a very economical solution, 
especially for strongly polarized beams, when we may get large rejection rates (roughly 
equal to the ratio of polarized to unpolarized Born cross section). The radical solution of 
this problem is to introduce longitudinal polarizations in the Born cross section as used 
in the crude and primary integrated cross section. At the moment they are completely 
unpolarized. 

In order to save CPU time a special method of "recycling" the r decay events is devised 
(see below). The interface supplies also the subroutine TRAL04, which is required by 
TAUOLA in order to transform r decay products to the laboratory frame. 

Some additional subroutines, necessary for the proper functioning of TAUOLA and 
PHOTOS, are placed in KK-all/KK2f/Tauface.f. 

• TaupairJnitialize(xpar) initializes TAUOLA and PHOTOS packages with the help of 
INIMAS, INITDK, INIPHY and DEKAY. It initializes the book-keeping for the spin 
weight, that is the weight used to introduce all spin effects in r decays. 

• Taupair_Finalize prints the average spin weight. 

• Taupair_Makel generates in the first step the unpolarized r decays using DEKAY. 
The polarimeter vectors m_HvecTaul and m_HvecTau2 are determined. 

• TaupairJmprintSpin introduces spin effects with the help of rejection using spin 
weight. The polarimeter vectors are sent to GPS with the help of GPS_SetHvectors 
and the spin weight is calculated with the help of GPS_MakeRho2; the event is then 
rejected or accepted. For the rejected event the r-pair event is "recycled", that 
is each r decay product is Euler-rotated and reused in the rejection method. The 
procedure is repeated until the event is accepted. The whole procedure is correct 
because we know exactly the average of the spin weight. 

• Taupair_Make2 transforms accepted r decay products to the CMS by calling DEKAY(ll) 
and DEKAY(12). The transformation is defined according to the GPS rules, for each 

t by GPS_tralorPrepare, and is performed with the help of GPS_TralorDolt hidden 
inside the TRAL04 routine. 

• Taupair_Clone performs the "recycling" of a r-pair by means of the Euler rotation 
in the rest frame of each r. 
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• Tralo4(Kto,P,Q,AM) , see above. 

• FILHEP(N,IST,ID,JM01,JM02,JDA1,JDA2,P4,PINV,PHFLAG) writes single parti- 
cles in r-decay into HepEvt class. For historical reasons HepEvt-Fill is not used 
directly. 

• Taupair_SetKeyClone(KeyClone) sets KeyClone. KeyClone switches between two op- 
erational modes of TaupairXlone. Both of the modes implement a valid solution. 

• Taupair_Getlslnitialized(lslnitialized) gets to know the outside world if TAUOLA is 
active (Islnitialized=l). 

• Taupair_GetHvectors(HvecFerl,HvecFer2) provides the polarimetric /i-vectors. 



4.13 Electroweak library 

The library of electroweak (EW) corrections is placed in the KK-all/dizet subdirectory. 
In the initialization phase the EW form factors dependent on s', fermion type and some of 
them (electroweak boxes) also on the scattering angles, are placed in the look-up tables. 
During the event generation they are interpolated in s' and cos 6 and provided to the 
GPS module, where the CEEX spin amplitudes are calculated, or used in the BornV_Dizet 
being used (through BornV_Differential) in the QED3, where the EEX distributions are 
calculated. The main aim of the above organization it to speed up the MC by using the 
EW form factor from look-up tables instead of calculating them, in fact many times, for 
each MC event. This reason may be even more important in the future version of the 
EW corrections which will be slower, owing to the inclusion of more genuine two-loop 
corrections. 

There are two modes of the initialization of the EW look-up tables. In the default 
mode they are calculated and stored in several disk files, each for one fermion type, and 
in the initialization of the /C/C MC run these tables are read by the BornV module. 
This mode is more conservative (safer), because the Fortran77 program providing elec- 
troweak corrections does not need to be linked and executed together with the proper 
MC event generator. We do not therefore need to worry about the clashes of the names 
of the procedures and common blocks, and the possible problems with re-initialization 
of the EW library for different types of the fermions is avoided^. In the actual im- 
plementation the EW library DIZET is run separately, under a special main program 
KK-all/dizet /TabMain. f, together with an interface module DZface, which acts as an 
interface to the Dizet library, properly setting the input data to Dizet and writing the 
EW form factors in the disk file. 

One important disadvantage of the above method is that the input parameters of the 
EW corrections, such as the Higgs mass, cannot be changed easily (for fitting), because 
it requires re-producing new look-up tables of the EW corrections. This is why we also 
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This reason seems to be now less important than in the early stages of the development of the /C/C 



MC, because most of the /C/C MC code now fulfils the rules of programming in sections 4.2 
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implemented the second interface to EW library DIZET in which the look-up tables are 
calculated in the initialization phase of the MC run (in this case the EW library is linked 
with the entire MC program). 

How are the EW tables produced in the default method? This is done by invoking in 
the KK-all/ dizet/ one of the commands: 

make tables # it makes all tables 

make table. mu # it makes . /table. mu using . /input. mu 

make table. tau # it makes . /table. tau using . /input. tau 

make table. down # it makes . /table . down using ./input. down 

make table. up # it makes ./table. up using ./input. up 

make table. botom # it makes . /table .botom using . /input .botom 

The input data should be used the same as in the MC run. In fact, the BornV module 
is checking if the important EW input data used to generate tables match the actual data 
provided by the user for the MC run. If not then the program stops. 

How does one avoid producing EW tables on the disk (and apply the second method)? 
For the instructions see KK-all/ff bench/Makefile. 

The interface DZface to DIZET of the Dubna-Zeuthen EWRC group version 6.21 
is based on the analogous interface in KORALZ 4.x. Let us now list and explain the 
subprograms in the interface module DZface 

• DZf ace_Initialize( KFfin, xpar) is the class initializer. Initialization of DIZET 
is done with CALL DIZET (NPAR, . . . ). NPAR and other input parameters are defined 
in the xpar vector, see Table B. 

• DZf ace_ReaDataX(DiskFile , iReset , imax, xpar) is functionally the same as 
KK2f_ReaDataX. 

• DZf ace_Tabluj fills in the EW form factors into look-up tables in /c_BornV/. The 
factors are provided by DZf ace_MakeGSW (see below). Tabulation is done in the three 
ranges of y/s with different numbers of points. The energy ranges and numbers of 
points are defined in BornV. h. 

• DZf ace_WriteFile(DiskFile) writes tables of the EW form factors into a disk file. 

• DZf ace_Clone (KFfin) copies tables calculated for the actual m_KFfin into tables 
for KFfin. It is used to create tables for c and s quarks using tables of u and d 
quarks (saving a little bit of CPU time). 

• DZf ace JVIakeGSW (Mode, ww,cosi,GSW, QCDcorN) gets the EW form-factors GSW and 
QCD corrections QCDcorN out of DIZET, at \fs =ww and cos# =cosi. EW form- 
factors are obtained with the CALL rokanc (...), while QCD corrections come from 
DZf ace_QCDtab (see below). 

• DZf ace_QCDtab (Mode, ww, QCDcorN) makes QCD corrections with the 

CALL qcdcof (...). It is done in an iterative way in order to find out the QCD 
corrections at a given \fs. 
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Note that the QED coupling constant the alf inv is separate from the alfinv used in 
the bremsstrahlung part of /C/C and the fermion masses in Dizet are isolated from those 
in /C/C. 

4.14 Random number generators 

The /C/C program in the present version uses exclusively the RAN MAR random num- 
ber generator P7|, jSSfl. It is reprogrammed as the pseudo-class PseuMar. Its single- 
precision generator is accessible with the help of the double-precision interface routine 
PseuMar_MakeVec. In this way we avoid possible interference with libraries of JETSET, 
PHOTOS and TAUOLA, which have their own independent random-number generators, 
and typically also have their own version of the RAN MAR generator. 

4.15 Other modules 

Every MC program of this size has to have its own tools for book-keeping of the MC 
weights and for making histograms of the weight distribution. In the /C/C MC the built-in 
histogramming package GLK plays this role. Histogramming entries are similar to those 
of the CERN library HBOOK. Apart from histogramming, it also has the capability of "mea- 
suring" several properties of the MC weight. This weight monitoring is done with the en- 
tries GLK_Mbook, GLK_Mf ill, GLK_Mget All , GLK_MgetNtot , GLK_MgetAve and GLK_Mprint. 
The GLK module features also simple, though versatile, graphical capabilities - it can 
plot histograms by exploiting the BTgXenvironment picture. 

The other auxiliary package is MathLib.f, which includes subprograms for Gaussian 
integration and some transcendental functions. The /C/C MC does not need any external 
mathematical library. 



5 Semianalytical "normalization anchor" 

In the complicated MC program aspiring to a precission of order 10 -4 , it is critical to 
master the overall normalization at this precision level. This can be done by comparing 
the program results with another MC program or with a semi-analytical calculation, that 
is the calculation in which as many integrations as possible are carried out analytically and 
the remaining ones are done with the Gaussian method^. Such a test of MC normalization 
is of critical importance - this is why we call it a normalization anchor. 

It is not trivial to analytically integrate the multiphoton phase space; nevertheless, 
for simple or simplified exponentiated distributions such as EEX and simple or absent 
kinematical cuts, such an integration is possible, see refs. [|l|,|39| . 



Here, since our aim is essentialy limited to the very precise numerical test of the 
MC phase-space integration, we have chosen the O(a ) EEX model, that is the Born 
differential cross section multiplied by the soft factors, see eqs. (|T5|) and ([[61). 



24 



This is the way the Zf itter phase space integration is done. 
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Figure 2: Test of the normalization of the low-level part of the Monte Carlo, for simplified QED 
multiphoton distribution. The difference of the MC result and semianalytical result divided by semi- 
analytical is plotted as a function of v maiX = 1 — s' min /s. Test is done for fi + pT final state at 
^fs = 189GeV. In case (a) t> max = 1 — 4m^/s is taken; the last bin represent the entire phase phase 
space. In case (b) v ma , x = 0.999. 



In this case it is relatively easy to obtain by analytical integration the O(a ) EEX 
cross section, keeping terms L Q a l , L 1 ^ 1 , L 2 a 2 , where L is the big logarithm, for both ISR 
and FSR. For the ISR, however, due to the Z resonance, and 7* peak at s' = 0, in order 
to reach the necessary precision of order 10~ 4 , it is necessary to analytically calculate two 
additional terms of orders L l a 2 and L 3 a 3 . This was done in ref. |12j| , and the relevant 
semianalytical formula looks as follows: 

f»mai 

°SAN = / dv ^BorntX 1 ~ ~ V )) D ISr(v) D FSR (u), 

JO 

— / 2\ / 1 

D FSR (u) = e^ 1 ^ e3 7/- 5 7/in(i-«)+f U+V; T/W t/- 1 M — ^ 7/ ln(l 



= e^+^ e * 7e+ -V2+^J 7e w 7e - 1 1 - -j e \n(l-v) - -- ln 2 (l - v) + 7 e 2 

V 4 2 7r 



(124) 



Note that the coefficient in front of the 0(L 3 a 3 ) term is zero, as marked explicitly. We 
have checked numerically that the ISR term of OiL^o?^) is worth several per cent for the 
cross section located close to v — 1. 

In fig. |2] we present the comparison of the JCJC MC with the above semianalytical 
formula. The difference between the MC result and the semianalytical result is divided by 
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the semianalytical result. The comparison is done for the final state at y/s =189 

GeV, as a function of v max . In the last point (bin) the entire phase space is covered, 
Vmax — 1 — 4m^/s. We conclude that we control the phase-space integration at the level 
of 2 x 1CT 4 for f max < 0.999, including the Z radiative return, and at the level of 3 x 1CT 3 
for no cuts at all. 



Parameter 



Position and meaning 



CMSene 

DelEne 

Ninp 

Nout 

LevPri 

IelPri 

Ie2Pri 

IdYFS* 

WtMax* 

KeyWgt 

IdeWgt* 

KeyELW 

vvmin* 
vvmax 

DelFac* 

NphMax** 

KeylSR 

KeyFSR 

KeyPia** 

mltlSR** 
mltFSR** 
KeyFix 

KeyWtm** 
KeyINT 

KeyGPS 

KeyQSR 



xpar(l) (=100): y/s, centre-of-mass (CMS) energy [GeV] 
xpar(2) (=0d0): Beam energy spread [GeV] 
xpar(3) (=5): Input unit number (unused) 
(=16): Output unit number 



xpar (4) 
xpar (5) 
xpar (6) 
xpar (7) 
xpar (8) 
xpar (9) 



=1 weight events 



(=0): PrintOut Level 0,1,2 
(=1): PrintOut Start point 
(=1): PrintOut End point 
(=600) : Pointer for internal histograms 
(=1): Maximum weight for rejection 
xpar (10) (=0): Switch between constant =0 and variable 
xpar (11) (=74): Ident of the EEX principal weight 

xpar (12) (=1) : Type of electroweak corrections, =0 only for tests, =1 default 
for DIZET 

xpar (16) (=lD-5): Minimum real photon energy in units of beam energy 
xpar (17) (=ld0): Maximum value of v = 1 — s'/s- variable, where s' is mass 
squared of // system including FSR photons! See more comments in the text, 
xpar (18) (=ld-3): FSR cut eps=vvmin*DelFac 

(=100): Hard-wired maximum photon multiplicity 



xpar (19) 
xpar (20) 
xpar(21) 
xpar (22) 



(=1) 
(=1) 
(=1) 



Test switch, KeyISR=0 swithes off the ISR 

Test switch, KeyFSR=0 switches off the FSR 

Removal of FSR photons below Emin=Ene*Delta in CMS, 



for KeyPia=0,l removal is OFF, ON 



xpar (23) (=0): 
xpar (24) (=0): 
xpar (25) (=0): 



: Special tests: fixed ISR multiplicity for mltISR>0 
: Special tests: fixed FSR multiplicity for mltFSR>0 
: Type of ISR, for KeyFix=0,l QED without beamstrahlung, 
for KeyFix=2 beamstrahlung is ON, see also KeyGrid 

xpar (26) (=0) : Special tests only: mass terms in "crude" MC photon distrib. 
xpar (27) (=2): Switch of ISR-FSR Interference (IFI), for KeyINT=0 it is 
OFF, for KeyINT=2 it is ON, KeyINT =1 is only for special tests 
xpar (28) (=1): Level of new exponentiation CEEX, note vmaxGPS overrules 
KeyGPS for each type of final fermion 

xpar (29) (=1): Photon emission from the final quarks is ON, OFF for 
KeyQSR=0,l 



Table 2: List of input parameters of the /C/C generator. General and related to QED radiation 
input parameters. Default values in brackets. User may change, with precautions, the starred 
items, while the doubly starred ones should never be changed. 
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Parameter 


Position and meaning 


Alflnv* 


xpar(30) (=137 . 0359895D0) : Inverse of oqed 


GNanoB* 


xpar(31) (=389.37966D3): GeV~ 2 — >[nb] conversion factor 


Gf ermi* 


xpar(32) (= 1.16639d-5): Fermi coupling [GeV~ 2 ] 


Technical cuts for non-IR real photon corrections etc. 


Xenph** 


xpar(40) (=1.25d0): Enhancement factor for "crude" photon multiplicity 


vcutl* 


xpar(41) (=l.d-9): Techn. cut for single non-IR bremss. correction 


vcut2* 


xpar(42) (=5.d-2): Techn. cut for double non-IR bremss. correction 


vcut3* 


xpar(43) (=0.1d0): Techn. cut for triple non-IR bremss. correction 


QCD flags/data 


KeyHad 


xpar (50) (=1) : Hadronization/showering flag, for KeyHad=0 , 1 hadronization 




is OFF,ON. 


HadMin* 


xpar (51) (=0.200d0): Minimum mass [GeV] for hadronization/showering 


alfQCD* 


xpar (52) (=0.118d0): a QC D 


Non zero beam polarization may require adjustment of WtMax 


spinlx 


xpar (61) (=0d0): polarization vector beam 1, x-component 


spinly 


xpar (62) (=0d0): polarization vector beam 1, y-component 


spinlz 


xpar (63) (=0d0): polarization vector beam 1, z-component 


spinlx 


xpar (64) (=0d0): polarization vector beam 2, x-component 


spinly 


xpar (65) (=0d0): polarization vector beam 2, y-component 


spinlz 


xpar (66) (=0d0): polarization vector beam 2, z-component 


Beamstrahlung parameters for Thorsten Ohl's package CIRCE 


IRCroots 


xpar (71) (=350d0): yfii [GeV] discrete values 350,500,800GeV 


IRCacc 


xpar (72) (=3d0): accelerator type 


IRCver 


xpar (73) (=5d0): version number 


IRCdat 


xpar (74) (=19980505d0) : date 


IRCxchat 


xpar (75) (=1): printout level 


KeyGrid* 


xpar (76) (=2): Activated by setting KeyFix=2 KeyGrid=2 invokes Foam, 




KeyGrid=0, — 1, +1 invokes Vegas: KeyGrid=0 creates and writes Vegas grid 




on the disk, KeyGrid=— 1 creates and dumps grid on the disk, =+1 reads only 



Table 3: List of input parameters of the /C/C generator in xpar vector. General parameters and 
QED radiation. Default values in brackets. User may change, with precautions, the starred 
items, while the doubly starred ones should never be changed. 



6 Use of the program 

In this section we will familiarize the reader with the input and output parameters, and 
the use of the present version of the /C/C Monte Carlo. We will present two simple 
demonstration main programs. Their double role is to serve as a useful template for the 
user to create his/her own main program and to help the user to check quickly that the 
newly installed /C/C generator runs correctly. We shall describe in detail all the input 
parameters of /C/C. 
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Parameter 


Position and meaning 


KFini* 


xpar(400) 


(= 11): Beam flavour code 




j'-th fermion is 


included in MC generation if its Mask(j)=l 


Mask( 1) 


xpar(401) 


(=1) 


Mask variable for d quark 


Mask( 2) 


xpar(402) 


( = 1) 


Mask variable for u quark 


Mask( 3) 


xpar(403) 


( = 1) 


Mask variable for s quark 


Mask( 4) 


xpar(404) 


( = 1) 


Mask variable for c quark 


Mask( 5) 


xpar(405) 


(=1) 


Mask variable for b quark 


Mask(13) 


xpar(413) 


(=1) 


Mask variable for muon lepton 


Mask(15) 


xpar415() 


(=1) 


Mask variable for tau lepton 


Basic electroweak input data 


MZ 


xpar(502) 


(=91 


. 187D0): Mass of Z-boson [GeV] (PDG 1996) 


SwSq 


xpar (503) 


(=.22276773D0): sm 2 (9 w ) where &w is EW mixing angle 


GammZ 


xpar(504) 


(= 2 


.50072032D0): Z width (from Dizet) 


MasPhot* 


xpar(510) 


(= 1D-60): Photon mass used as IR regulator 


The data base record below is for d quark, j = 1 


KFferm(j)* 


xpar(501+10*j) 


(= 1): Flavour code 


NCf (j)* 


xpar(502+10*j) 


(= 3) : Number of colours 


Qf (j)* 


xpar(503+10*j) 


(=-1): 3 x charge 


T3f (j)* 


xpar(504+10*j) 


(=-1): 2xT3L =2xIsospin for left component 


Helic(j)* 


xpar(505+10*j) 


(= 1): 2xhelicity, not used 


Mferm(j)* 


xpar(506+10*j) 


(= O.OlOdO): Mass [GeV] (PDG) 


MfCon(j)* 


xpar(506+10*j) 


(= O.lOOdO): Constituent mass, not used 


WtMax(j)* 


xpar(507+10*j) 


(= 5.0d0): Maximum weight for rejection 


AuxPar(j)* 


xpar(508+10*j) 


(= 0.99d0): below vmaxGPS CEEX, above EEX 


More electroweak input data 


Ibox 


xpar (801) 


(= 1): EW box flag, input for Dizet 


MH 


xpar (805) 


(= 100D0): Higgs mass, input for Dizet 


Mtop 


xpar(806) 


(= 175D0): Top mass, input for Dizet 


NPAR(j) 


xpar ( 900+ j) ,j = 


1,21: Input flags for Dizet 6.21 



Table 4: List of input parameters of the /C/C generator. Initial/final fermion properties and 
EW parameters. Default values in brackets. User may change, with precautions, the starred 
items, while the doubly starred ones should never be changed. 



6.1 Principal entries of /C/C 

The principal entries of the /C/C package, which the user will call in his/her application 
in order to generate a series of MC events, were already listed and described briefly in 
Section [|. Here we shall add more information on their functionality. The calling sequence 
constituting a typical Monte Carlo run will look as follows: 

CALL KK2f_ReaDataX(' . / .KK2f \_def aults ' ,1,10000 , xpar) ! reading default input 

CALL KK2f_ReaDataXC ./user. input' ,0, 10000, xpar) ! reading user's input 

CALL KK2f _Initialize (xpar) ! initialize generator 

DO loop=l, 10000 ! loop over MC events 
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Parameter 



Position and meaning 



Jakl 

Jak2 

idff** 

itdkRC 

xkOdec* 

KeyAl 

Cabib** 

GV* 

GA* 



xpar(2001) (=0): First r decay mask 

xpar(2002) (=0): Second r decay mask 

xpar(2003) (=15): PDG ident of the first r 

xpar(2004) (=1): QED rad. switch in leptonic decays 

xpar(2005) (=0.001d0): IR-cut for QED rad. in leptonic decays 

xpar(2006) (=ld0): Type of a\ current 

xpar(2007) (=0.975d0): Cosine of Cabibbo angle 

xpar(2008) (= IdO): Vector coupling gy in r decay 

xpar(2009) (=-ld0): Axial coupling g& in r decay 



xpar(2010) (= 0.5d0): In 3-pion decay BR of 7r + 7r"7r~ (vs 7r~7r°7r°) 

xpar(2011) (=0.6667d0): In K* decay BR of K + 7t° (vs tt + K°) 

xpar(2012) (=0.5d0): Probability of K° to be K s 

xpar(2013) (=0.5d0): Probability of K% to be K s 



BRA1 
BRKS 
BRKO 
BRKOB 



Branching ratios 



BRAE 



BRAMU 
BRAPI 
BRA2PI 
BRA3PI 



xpar(2101) (=17 . 810d-2) : Branching ratio r~ e~ . IMPORTANT! Entry 

2101 set smaller than -IdO will activate internal defaults of Tauola. In such a 

case all input from 2008-2122 will be IGNORED 

xpar(2102) (=17 . 370d-2) : Branching ratio r~ — > \T 

xpar(2103) (=11 . 080d-2) : Branching ratio r~ — > 7r~ 

xpar(2104) (=25 . 320d-2) : Branching ratio r~ — > 7r~,7r° 

xpar(2105) (=18 . 380d-2) : Branching ratio r~ — > a{ 



Other branching ratios are in xpar (2106-2122), see .KK2f_defaults 



Table 5: Input parameters for the TAUOLA package. For a complete description, see the 
Manual of TAUOLA Default values in brackets. User may change, with precautions, the 
starred items, while the doubly starred ones should never be changed. 



CALL KK2f_Make ! generate single event 
ENDD0 

CALL KK2f _Finalize ! final book-keeping 

CALL KK2f _GetXSecMC (XSecPb , XErrPb) ! get total cross section 

In the first call of KK2f _ReaDataX, default data are read into the array REAL*8 xpar (10000) 
The /C/C itself has almost no data hidden in the source code. (This is not true for TAUOLA 
and JETSET). The file . KK2f _def aults is read first into array xpar. This file of defaults 
is provided in the main distribution directory. The user should never modify it. It can 
be copied to a local directory or, better, a symbolic link should be created to the orig- 
inal file. The . KK2f .defaults is rather large and the user is usually interested only in 
changing some subset of these data. In the second call on KK2f _ReaDataX the user can 
overwrite the default data with his/her own smaller set of input data, which are placed 
in the user . input file. See next subsection for more details on the input data. 

The KK2f .Initialize is invoked to initialize the generator. It reads input data from 
array xpar, prints them and sends them down to various modules and auxiliary libraries. 
The program entries have to be called in strictly the same order as in the above example. 
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At this point we are ready to generate a series of the MC events. The generation of a 
single event is done with the help of KK2f _Make. After the generation loop is completed, 
we may invoke KK2f .Finalize, which does final book-keeping, prints various pieces of 
information on the MC run, and calculates the total MC integrated cross section and its 
statistical error in units of picobarn. In order to access the total cross section the user 
should call the routine KK2f_GetXSecMC(XSecMC,XErrMC). 



6.2 Input data 

As we stated previously, in the second call on KK2f _ReaDataX the users can overwrite the 
default with their own preferred values. Note that the user should never modify certain 
data items (without consulting authors of the program) and that the other ones can be 
changed by the user, see below. For example the simplest input data, which define only 
the CMS energy, look as follows: 

BeginX 

*<ia>< data >< comments > 

1 190d0 CmsEne =CMS total energy [GeV] 

EndX 

As we see, data cards start with the keyword BeginX and end with the keyword EndX. The 
comment lines are allowed - they start with * in the first column. In the comments we 
specify the meaning of the data, their range, and whether the user is allowed to modify 
them. The data themselves are in a fixed format, with the address i in xpar(i) followed 
by the data value and trailing comment. The four examples of input data sets for the two 
demonstration programs f f bench/demo . f and f f bench/ProdMC . f in the subdirectories 
ffbench/Mu, ff bench/Inclusive and the other ones, provide useful templates for the 
typical user data. The complete set of all user data in KK2f _def aults is described in 
detail in Tables Understandably, the user will manipulate, in most cases, only a 

small subset of the data and, in most cases, will stick to the default values. 



6.3 MC events and other output 



The principal output of /C/C is the Monte Carlo event, which is just a list of final-state 
four-momenta in [GeV] units and flavours, encoded in the standard /d_HepEvt/ common 
block, see section |4~4[ 

All beam, photon and parton momenta before hadronization are available alternatively 



through "getter" subroutines from class KK2f, see section 4.3 or HepEvt: 



DOUBLE PRECISION 
INTEGER 



pi (4) ,p2(4) ,p3(4) ,p4(4) ,PhoAll(100,4) 
NphAll 



CALL HepEvt_GetBeams(pl ,p2) 
CALL HepEvt_GetFf ins(p3,p4) 
CALL HepEvt_GetPhot All (NphAll , PhoAll) 



get beam momenta 

get momenta of two final fermions 
get photon multiplicity and momenta 
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where NphAll is the total photon multiplicity (see also the f f bench/ProdMC . f example). 
Note that beamstrahlung photons are added to the record as two zero-angle ISR photons, 
so that total energy is conserved. Alternatively, beamstrahlung photon momenta are also 
available through a dedicated getter: 

DOUBLE PRECISION PhoBst(100,4) 
CALL HepEvt_GetPhotBst (NphBst , PhoBst) 



Parameter 



Position and meaning 



EEX C(a u ) 
EEX Oia 1 ) 
EEX C(q 2 ) 
EEX <D(a 3 ) 
CEEX 0(a°) 
CEEX Oia 1 ) 
CEEX 0(a 2 ) 
CEEX O(a ) 
CEEX dice 1 ) 
CEEX 0{a 2 ) 



WtSet(71) 

WtSet(72) 

WtSet(73) 

WtSet(74) 

WtSet(201) 

WtSet(202) 

WtSet(203) 

WtSet(251) 

WtSet(252) 

WtSet(253) 



without ISR-FSR interference 
without ISR-FSR interference 
without ISR-FSR interference 



Table 6: The meaning of the weights in the WtSet. 



6.4 Weighted events, alternative weights 

Normally, the user will run the program in the mode with the weight equal to 1. Running 
in the mode with weighted events may be useful for various tests. It can be useful, for 
example as a cross check, in the situation when one selects output events strongly, that 
is imposes cuts that eliminate all but say 1 event in a 1000. If at the same time it is 
seen from the output of KK2f .finalize that the cross section corresponding to w > u> max 
is at the similar level of 10~ 3 , it is then necessary to cross check if the accepted events 
do not coincide, by bad luck, with the "overweighted" events. If it were true, then the 
cross section and the distribution of the accepted events could be affected by factor of 2 
or more. In that sense the weighted events are "safer". 

It should be kept in mind that, although we have set the maximum weights for the 
rejection rather high, the user may try an untested configuration of the input data for 
which the cross section corresponding to w > u> max is too high. We recommend that the 
user always check, at the end of the run, the output from KK2f _f inalize, the table in 
which the percentage of the "spill over" cross section corresponding to w > w max is given. 

The other advantage of the weighted events is that in most cases one needs less CPU 
time to get the same statistical error in the cross sections and in the histogram. It can 
be profitable if one needs to perform many runs with various input parameters. 

In the run with weighted events the user may access the principal weight WtMain and 
the auxiliary weights WtSet through another getter: 
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DOUBLE PRECISION WtSet(lOOO), WtMain, WtCrud 
CALL KK2f _GetWt All (WtMain , WtCrud , WtSet ) 

see also the f f bench/ProdMC . f example. The actual auxiliary weight should be defined as 
WtCrud*WtSet (i) . Note that events with WtCrud=0d0 may have undefined four- momenta, 
so the user should protect his program against crashing upon an attempt of working out 
the kinematics of such an event. 

The weights WtCrud and WtSet are also defined in the run with WtMain=ldO. They 
can be recorded and used in the subsequent run in order to estimate the effects that 
are included or excluded in the auxiliary weight WtSet (i). The meaning of the most 
important weights in WtSet is described in table |6|. 

How to get cross sections and distributions corresponding to WtSet (i) using an event 
generated with WtMain =1 and recorded on the tape? 

• The user should record on the tape the vector WtSet for each event, together with 
WtCrud and WtMain. 

• In the subsequent run, to weight events from the tape, each event should be weighted 
with the ratio0 WtSet (i) / WtSet (203), because in the standard case 
WtMain=WtCrud*WtSet (203) . 

A typical application of the above method could be to find out, for a given arbitrary 
distribution or cross section, the estimate of physical precision due to higher orders. We 
recommend that the user take half of the difference 0(a 2 ) — 0{a 1 ) as an estimate of 
the physical precision. This can be calculated by applying the above method with the 
following weight: (WtSet (203) -WtSet (202) ) *WtCrud/Wtmain. 

The above method cannot be used for varying the input parameters of the SM, such 
as the Higgs mass, because this would require recalculating WtSet. We may provide such 
a capability in the next versions. 

6.5 Frequently asked questions on program use 

Some additional information, useful for practical use of the program is collected as answers 
to "frequently asked questions": 

• How does one properly normalize total cross section? Look into two demonstration 
programs in ff bench subdirectory. 

• How does one update tables of electroweak corrections? 

cd KK-all/dizet 

make all (or make table. tau etc.) 

25 This will not work for muons with v > 0.999 and quarks with v > 0.99 where, for technical reasons 
we use WtMain=WtCrud*WtSet(74) . This restriction is not important for most of practical purposes. 
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For more details see section |4.13| . 

• How does one switch on beamstrahlung? Include KeyFix=2 and KeyGrid=2 in the 
user input data. An example program is included in ff bench subdirectory. 

• How does one switch off radiation for quarks ? Include KeyQSR=0 in the user input 
data. 

• How does one switch from CEEX to EEX for quarks ? Include vmaxGPS=0 for all 
quarks in the user input data. 

• How does one update compilation flags everywhere? Compilation flags are set for 
AIX. Examples of f77 flags for HPUX, Linux, ALPHA are in . /ff bench/Makefile. 
In order to update centraly makefiles in all subdirectories do the following: 

cd ./ff bench 
make makflag 

This causes the mapping makefile .template — > makefile in all subdirectories, 
updating compilation flags everywhere with the ones from ./ffbench/Makefile. 

• How does one calculate the QED physical error for a given observable? Calculate the 
difference between C(« 2 )ceex and C(q; 1 )ceex and the difference between C(a 2 )cEEX 
and 0(« 3 )eex- This can be done by running the MC with weighted events and 
taking the difference of the weights or with unweighted events, following instructions 
in the previous section. 

• How can one be sure about the technical precision? The problem may arise for strong 
selection cuts. In this case we advise the user to rerun the program with weighted 
events and check whether the results are the same. 



7 Outlook and conclusions 

As is summarized in Table [I], the present version of the /C/C MC has almost the full 
functionality of the older KORALZ and KORALB event generators. The most important 
new features in the present /C/C are the ISR-FSR interference, the second-order subleading 
corrections, and the exact matrix element for two hard photons. This makes /C/C already 
a unique source of SM predictions for the LEP2 physics program. The inclusion of the 
beamstrahlung makes it useful for the LC studies. Note that for these the electroweak 
correction library has to be reexamined. The most important omission in the present 
version is the lack of neutrino and electron channels. Let us stress that the present 
program is an excellent starting platform for the construction of the second-order Bhabha 
MC generator based on CEEX exponentiation. We hope to be able to include the Bhabha 
and neutrino channels soon, possibly in the next version. The other important directions 
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for the development are the inclusion of the exact matrix element for three hard photons, 
together with virtual corrections up to 0(a 3 L 3 ) and the emission of the light fermion 
pairs. The inclusion of the W + W~ and tt final states is still in a farther perspective. 
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Appendix: Output of the demonstration program 



Demo for KK MC 



1000 requested events 

*******************************************:+;** 



sic*** **** 
********** 
******* 
********** 
**** ***** 
**** ***** 



**** **** 
**** **** 
**** **** 
********* 
****** 
******** 
**** **** 
**** **** 



*** *** 
**** **** 
***** ***** 
************* 
************* 
**** *** **** 
**** * **** 
**** **** 



********************************************** 
********************************************************* 

* KK Monte Carlo 

Version 4.13 25 Jan. 2000 



************** 
****** 
********** 
***** *** 

***** *** 
********** 
******* 
************** 
************** 



* 


200.00000000 


CMS energy average 


CMSene 


al 


* 


* 


.00000000 


Beam energy spread 


DelEne 


a2 


* 


* 


100 


Max. photon mult. 


npmax 


a3 


* 


* 


1 


ISR switch 


KeylSR 


a4 


* 


* 


1 


FSR switch 


KeyFSR 


a5 


* 


* 


2 


ISR/FSR interferenc 


KeyINT 


a6 


* 


* 


1 


New exponentiation 


KeyGPS 


a7 


* 


* 


1 


Hadroniz. switch 


KeyHad 


a7 


* 


* 


.20000000 


Hadroniz. min. mass 


HadMin 


a9 


* 


* 


1.00000000 


Maximum weight 


WTmax 


alO 


* 


* 


100 


Max. photon mult. 


npmax 


all 


* 


* 


11 


Beam ident 


KFini 


al2 


* 


* 


.00100000 


Manimum phot. ener. 


Ene 


al3 


* 


* 


. 10000000E-59 


Phot. mass, IR regul 


MasPho 


al4 


* 


* 


1.2500000 


Phot . mult . enhanc . 


Xenph 


al5 


* 


* 


. 10000000E-08 


Vcutl 


Vcutl 


al6 


* 


* 


.50000000E-01 


Vcut2 


Vcut2 


al6 


* 


* 


.00000000E+00 


Vcut3 


Vcut2 


al6 


* 


* 


.00000000 


PolBeaml(l) 


Pollx 


al7 


* 


* 


.00000000 


PolBeaml(2) 


Polly 


al8 


* 


* 


.00000000 


PolBeaml(3) 


Pollz 


al9 


* 


* 


.00000000 


PolBeam2(l) 


Pol2x 


a20 


* 


* 


.00000000 


PolBeam2(2) 


Pol2y 


a21 


* 


* 


.00000000 


PolBeam2(3) 


Pol2z 


a22 


* 


*************************************************************************** 


*************************************************************************** 


* 


BornV Initializator 








* 


* 


91.18700000 


Z mass [GeV] 


amz 


al 


* 


* 


100.00000000 


Higgs mass [GeV] 


amh 


a2 


* 


* 


175.00000000 


Top mass [GeV] 


amtop 


a3 


* 


* 


2.50072032 


Z width [GeV] 


gammz 


a4 


* 


* 


.22276773 


sin(theta_w) **2 


sinw2 


a5 


* 


* 


137.03598950 


l/alfa_QED at Q=0 


Alf Inv 


a6 


* 


* 


.20000000 


MassCut light qqbar 


HadMin 


a6 


* 


* 


11 


KF code of beam 


KFini 


a7 


* 


* 


1.0000000 


Input vvmax 


vvmax 


a8 


* 


* 


.99999888 


reduced vvmax in MC 


vvmax 


a9 


* 


* 


Test switches: 








* 


* 


1 


Electroweak lib. 


KeyElw 


10 


* 


* 


1 


Z on/off switch 


KeyZet 


11 


* 


* 





mass terms on/off 


KeyWtm 


12 


* 


******************************************************************* 


*************************************************************************** 


* 


BornV Reading from disk 


file: 






* 


* 


. . / . . /dizet/table . down . 340pt 






* 


* 


91.18700000 


Z mass 


amz 


al 


* 


* 


100.00000000 


Higgs mass 


amh 


a2 


* 


* 


175.00000000 


Top mass 


amtop 


a3 


* 


* 


.22302485 


sin**2(thetaW) 


swsq 


a3 


* 


* 


2.49925439 


Z width 


gammz 


a3 


* 


* 


80.37787000 


W mass 


amw 


a3 


* 


* 


2.08825838 


W width 


gammw 


a3 


* 



**************************************************** 

**************************************************** 



****************** 



* 




BornV Reading from disk 


file: 






* 


* 




. . / . . /dizet/table . up . 340pt 






* 


* 


91 


. 18700000 


Z mass 


amz 


al 


* 


* 


100. 


,00000000 


Higgs mass 


amh 


a2 


* 


* 


175 


.00000000 


Top mass 


amtop 


a3 


* 


* 




.22302485 


sin**2(thetaW) 


swsq 


a3 


* 


* 


2. 


.49925439 


Z width 


gammz 


a3 


* 


* 


80 


,37787000 


W mass 


amw 


a3 


* 


* 


2. 


.08825838 


W width 


gammw 


a3 


* 



***** 

***** 



**************************************************** 

**************************************************** 
BornV Reading from disk file: 
. . / . . /dizet/table . down . 340pt 
91.18700000 Z mass 

100.00000000 Higgs mass 

175.00000000 Top mass 

.22302485 sin**2 (thetaW) 

2.49925439 Z width 



****************** 
****************** 







* 






* 


amz 


al 


* 


amh 


a2 


* 


amtop 


a3 


* 


swsq 


a3 


* 


gammz 


a3 


* 
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* 80.37787000 W mass amw a3 * 

* 2.08825838 W width gammw a3 * 
*************************************************************************** 
********************************************** 



* 




BornV Reading from disk file: 






* 


* 




. . / . . /dizet /table . up . 340pt 






* 


* 


91 


. 18700000 


Z mass 


amz 


al 


* 


* 


100. 


,00000000 


Higgs mass 


amh 


a2 


* 


* 


175 


.00000000 


Top mass 


amtop 


a3 


* 


* 




.22302485 


sin**2(thetaW) 


swsq 


a3 


* 


* 


2. 


,49925439 


Z width 


gammz 


a3 


* 


* 


80. 


,37787000 


W mass 


amw 


a3 


* 


* 


2. 


,08825838 


W width 


gammw 


a3 


* 



*************************************************************************** 

*************************************************************************** 



* 




BornV Reading 


from disk file: 






* 


* 




. . / . . /dizet/table . botom . 340pt 






* 


* 


91 


. 18700000 


Z mass 


amz 


al 


* 


* 


100. 


,00000000 


Higgs mass 


amh 


a2 


* 


* 


175. 


,00000000 


Top mass 


amtop 


a3 


* 


* 




,22302485 


sin**2(thetaW) 


swsq 


a3 


* 


* 


2. 


,49925439 


Z width 


gammz 


a3 


* 


* 


80 


,37787000 


W mass 


amw 


a3 


* 


* 


2. 


,08825838 


W width 


gammw 


a3 


* 



* 




BornV Reading from disk file: 






* 


* 




, . / . . /dizet/table 


. mu . 340pt 






* 


* 


91 


. 18700000 


Z mass 


amz 


al 


* 


* 


100. 


,00000000 


Higgs mass 


amh 


a2 


* 


* 


175. 


,00000000 


Top mass 


amtop 


a3 


* 


* 




,22302485 


sin**2(thetaW) 


swsq 


a3 


* 


* 


2. 


,49925439 


Z width 


gammz 


a3 


* 


* 


80 


,37787000 


W mass 


amw 


a3 


* 


* 


2. 


,08825838 


W width 


gammw 


a3 


* 



*************************************************************************** 
*************************************************************************** 



* 




BornV Reading 


from disk file: 






* 


* 




, . / . . /dizet/table . tau . 340pt 






* 


* 


91 


. 18700000 


Z mass 


amz 


al 


* 


* 


100. 


,00000000 


Higgs mass 


amh 


a2 


* 


* 


175. 


,00000000 


Top mass 


amtop 


a3 


* 


* 




,22302485 


sin**2(thetaW) 


swsq 


a3 


* 


* 


2. 


,49925439 


Z width 


gammz 


a3 


* 


* 


80 


,37787000 


W mass 


amw 


a3 


* 


* 


2. 


,08825838 


W width 


gammw 


a3 


* 



*************************************************************************** 

*************************************************************************** 



* 


KarLud_Initialize START 




* 


* 200.00000000 


CMS energy average 


CMSene 


== * 


* .00000000 


Beam energy spread 


DelEne 


== * 


* 1 


ISR on/off switch 


KeylSR 


== * 


* 


Type of ISR 


KeyFix 


== * 


* 1 


Elect_weak switch 


KeyZet 


== * 


* 


Fixed nphot mult . 


MltlSR 


== * 


* 50 


Max. photon mult. 


nmax 


== * 


*************************************************************************** 


* 4370.16701351 


xs_crude vesko 


xcvesk 


* 


* 4363.90727756 


xs_crude gauss 


xcgaus 


* 


* .00143443 


xcvesk/xcgaus-1 




* 




* 


KarLud_Initialize END 




* 



*************************************************************************** 



* 






KarFin 


Initialize START 






* 


* 




1 




FSR radiation on/off 


KeyFSR 


al 


* 


* 




1 




radiation from quark 


KeyQSR 


a2 


* 


* 




1 




removal switch 


KeyPia 


a3 


* 


* 


. 10000000E- 


-02 




infrared cut FACTOR 


delf ac 


a4 


* 


* 


. 10000000E- 


-07 




infrared cut itself 


delta 


a5 


* 


* 


. 10000000E- 


-02 




EminCMS for removal 


[GeV] 


a6 


* 


* 




50 




Max. photon mult. 


nmax 


a7 


* 


* 






KarFin 


Initialize END 






* 



*************************************************************************** 



* 




GPS Initializator 








* 


* 


91. 


, 18700000 


Z mass [GeV] 


MZ 


al 


* 


* 


2 


,49925439 


Z width [GeV] 


GammZ 


a2 


* 


* 




,22302485 


sin(theta w)**2 


Sw2 


a3 


* 


* 


137 i 


,03598950 


l/alfa_QED at Q=0 


Alf Inv 


a4 


* 


* 




Test switches: 








* 


* 




1 


Z on/off switch 


KeyZet 


a5 


* 


* 




1 


Electroweak lib. 


KeyElw 


a6 


* 


* 




1 


CEEX level 


KeyGPS 


a7 


* 


* 




1 


ISR emission 


KeylSR 


a8 


* 


* 




1 


FSR emission 


KeyFSR 


a9 


* 


* 




2 


ISR*FSR interferenc 


KeyINT 


alO 


* 
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* 






KK2f : 


Initialization 






* 


* 


9.4892226 






x-crude [nb] 


Xcrunb 


** 


* 


* 


List 


of final 


f ermions : 








* 


* 




1 




KF of final fermion 


KFf in 


** 


* 


* 


. 10000000 






mass of final ferm. 


amf erm 


** 


* 


* 


1.3410294 






Xborn [R] 


Xborn 


** 


* 


* 


5.0000000 






WtMax sampling par. 


WtMax 


** 


* 


* 


.99000000 






vmax for CEEX 


vmaxGPS 


** 


* 


* 




2 




KF of final fermion 


KFf in 


** 


* 


* 


. 10000000 






mass of final ferm. 


amf erm 


** 


* 


* 


2. 1445691 






Xborn [R] 


Xborn 


** 


* 


* 


5.0000000 






WtMax sampling par. 


WtMax 


** 


* 


* 


.99000000 






vmax for CEEX 


vmaxGPS 


** 


* 


* 




3 




KF of final fermion 


KFf in 


** 


* 


* 


. 20000000 






mass of final ferm. 


amf erm 


** 


* 


* 


1.3410294 






Xborn [R] 


Xborn 


** 


* 


* 


5.0000000 






WtMax sampling par. 


WtMax 


** 


* 


* 


.99000000 






vmax for CEEX 


vmaxGPS 


** 


* 


* 




4 




KF of final fermion 


KFf in 


** 


* 


* 


1.3000000 






mass of final ferm. 


amf erm 


** 


* 


* 


2.1445691 






Xborn [R] 


Xborn 


** 


* 


* 


5.0000000 






WtMax sampling par. 


WtMax 


** 


* 


* 


.99000000 






vmax for CEEX 


vmaxGPS 


** 


* 


* 




5 




KF of final fermion 


KFf in 


** 


* 


* 


4.5000000 






mass of final ferm. 


amf erm 


** 


* 


* 


1.3410294 






Xborn [R] 


Xborn 


** 


* 


* 


5.0000000 






WtMax sampling par. 


WtMax 


** 


* 


* 


.99000000 






vmax for CEEX 


vmaxGPS 


** 


* 


* 




13 




KF of final fermion 


KFf in 


** 


* 


* 


. 10565830 






mass of final ferm. 


amf erm 


** 


* 


* 


1.2225177 






Xborn [R] 


Xborn 


** 


* 


* 


8.0000000 






WtMax sampling par. 


WtMax 


** 


* 


* 


.99900000 






vmax for CEEX 


vmaxGPS 


** 


* 


* 




15 




KF of final fermion 


KFf in 


** 


* 


* 


1.7770000 






mass of final ferm. 


amf erm 


** 


* 


* 


1.2225177 






Xborn [R] 


Xborn 


** 


* 


* 


8.0000000 






WtMax sampling par. 


WtMax 


** 


* 


* 


1.0000000 






vmax for CEEX 


vmaxGPS 


** 


* 



******************************************************************** 



*************************************************************************** 

* KK interface of Tauola * 

* 2 Cloning procedure KeyClone tOl * 
*************************************************************************** 
*************************************************************************** 

* Parameters passed from KK to Tauola: * 

* dec. type 1-st tau Jakl tOl * 

* dec. type 2-nd tau Jak2 t02 * 

* 1 current type al dec. KeyAl t03 * 

* 15 PDG id 1-st tau idff t04 * 

* 1 R.c. switch lept dec itdkRC t05 * 

* . 10000000E-02 IR-cut for lept r.c. xkOdec t06 * 
*************************************************************************** 
*************************************************************************** 

* TAUOLA Initialization SUBROUTINE INIMAS: * 

* Adopted to read from KK * 

* 1.7770000 AMTAU tau-mass **** *** * 

* .51099900E-03 AMEL electron-mass **** *** * 

* . 10565830 AMMU muon-mass **** *** * 
*************************************************************************** 

skipped output from TAUOLA 

*************************************************************************** 

* *****TAU0LA LIBRARY: VERSION 2.6 ****** * 

* ***********August 1995*************** * 

* **AUTH0RS: S.JADACH, Z.WAS***********»* * 

* **R. DECKER, M. JEZABEK, J.H.KUEHN***** * 

* ** AVAILABLE FROM: WASM AT CERNVM ****** * 

* ***** PUBLISHED IN COMP. PHYS. COMM.*** * 

* *******CERN-TH-5856 SEPTEMBER 1990***** * 

* *******CERN-TH-6195 SEPTEMBER 1991***** * 

* *******CERN TH-6793 NOVEMBER 1992***** * 

* **5 or more pi dec: precision limited * 

* ****DEKAY ROUTINE: INITIALIZATION****** * 

* JAK1 DECAY MODE TAU+ * 

* JAK2 = DECAY MODE TAU- * 
*************************************************************************** 

*************************************************************************** 

Event listing (summary) 



I 


particle/jet KS 


KF 


orig 


P- 


X 


P- 


-y 


P- 


z 


E 


m 


1 


!e-! 


21 


11 







.000 




.000 


100 


.000 


100 


.000 


.001 


2 


!e+! 


21 


-11 







.000 




.000 


-100. 


.000 


100. 


,000 


.001 


3 


(tau-) 


11 


15 


1 


-11 


.368 


-8 


.313 


-16. 


.590 


21. 


,834 


1.777 


4 


(tau+) 


11 


-15 


1 


77 


.073 


19 


.642 


34. 


.404 


86. 


,677 


1.777 


5 


gamma 


1 


22 


1 


-65 


.705 


-11. 


,329 


6 


.665 


67. 


,007 


.000 


6 


gamma 
gamma 


1 


22 


1 




.001 




.000 


-24 


.480 


24. 


,480 


.000 


7 


1 


22 


2 




.002 




.001 




.000 




.002 


.000 


8 


nu tau 


1 


16 


3 




.030 




.121 




.003 




.125 


.010 


9 


(rho-) 


11 


-213 


3 


-11 


.338 


-8 


.192 


-16. 


.587 


21 


,709 


.690 


10 


pi- 


1 


-211 


9 


-6 


.338 


-4 


.212 


-9. 


.177 


11. 


,922 


.140 


11 


piO 


1 


111 


9 


-5 


.000 


-3 


.980 


-7 


.410 


9 


.787 


.135 


12 


nu tau" 


1 


-16 


4 


34 


.604 


8 


.167 


14. 


.712 


38 


,478 


.010 


13 


pi+ 


1 


211 


4 


42 


.469 


11 


.475 


19. 


.692 


48. 


,199 


.140 




sum: 


.00 






.000 




.000 




.000 


200. 


,000 


200.000 



* KarLud final report * 
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* 90875 

* 4370.16701351 

* 4363.65383283 

* 4370.1670 

* .00000000E+00 
******************* 
******************* 
* 

* 90875 

* 

* .01551546 

* 67.80514951 



total no of events nevtot 

ISRcru [R] ISRcru 

.28966439 ISRbest [R] , ISRerr ISRbest 

XKarlud [R] XKarlud 

KError [R] KError 

ijc***^*************************:^**********^****** 
***************************************** 

Report on wt_ISR of KarLud 

total no of events 
wt<0 events 
+- .00666845 <wt> 
+- .45215550 sigma of KarLud [R] 
*************************************************************************** 
*************************************************************************** 

* .99850963 +- .00006638 Average WT of Veskl AVeskl == * 

* 4363.90727756 +- .04363907 xs_est gauss [R] xcgaus == * 

* -.00005808 +- .00006738 xcve/xcgs-1 == * 
*************************************************************************** 
*************************************************************************** 



nevtot 
nevneg 
wt_ISR 
xskarl 



== * 
== * 



== * 
== * 



* 






KarFin Finalize START 






* 


* 


180070 




generated events 


nevgen 


a2 


* 


* 


.99710113 


+- 


.00017970 kinematics, smin 


wtl 


a5 


* 


* 


.99980533 


+- 


.00002691 j aerobian 


wt2 


a6 


* 


* 


.99929387 


+- 


.00061237 photon ang. dist . 


wt3 


a7 


* 


* 






ON MASS WEIGHTS 






* 


* 


.95525378 


+- 


.00052108 removal wgt wtrem 




bl 


* 


* 


89774 




no. of raw events 




b2 


* 


* 







wt6=0 events 




b3 


* 


* 


.99957099 


+- 


.00051944 control wgt wctrl 




b4 


* 


* 







marked photons 


MarTot 


a5 


* 


* 


. 10000000E-02 




emin 




b6 


* 


* 


. 10000000E-07 




delta 




b7 


* 


* 


. 15689398 




raw ph. multipl. 




b8 


* 


* 


6.00000000 




Highest phot. mult. 




b9 


* 


* 






YFSfin Finalize END 






* 



********************** 
********************** 



************************************* 
KarFin Finalize 



******** 



******** 



* 180070 

* .99634044 +- .00064027 

* .15643916 

*********************************************************** 



generated events 
general weight 
aver. ph. multi. 



-640 KK2f : Photon raw multiplicity 

nent sum bmin bmax 

1000 .00000E+00 .00000E+00 .33000E+03 

undf ovef sumw avex 

.00000E+00 .00000E+00 .21210E+04 .28322E+01 



nevgen 
wt 

avmlt 
**************** 



a2 
al 
a3 



0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 



.600000D+02 OXXXXXXXXXXXX I 
. 277000D+03 OXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX I 
.330000D+03 OXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 



. 196000D+03 OXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
. 990000D+02 OXXXXXXXXXXXXXXXXXXX 
.290000D+02 0XXXXX 
.800000D+01 OX 
.100000D+01 
.000000D+00 

,****2*2*$* ) J*2****» t ********** :*****: *„ ********** *********** M 

KK2f .Finalize printouts * 

cms energy total cmsene aO * 

total no of events nevgen al * 

principal info on x-section ** * 

+- .51653089 xs_tot MC R-units xsmc al * 

xs_tot picob. xSecPb a3 * 

error picob. xErrPb a4 * 

relative error erel a5 * 

WTsup, largest WT WTsup alO * 

** some auxiliary info ** * 

xs_born picobarns xborn all * 

Raw phot, multipl. === * 

Highest phot. mult. === * 

End of KK2f Finalize * 

*************************************************************************** 

skipping some lines 

*************************************************************************** 



* 






* 


200. 


.00000000 


* 




1000 


* 




** 


* 


47 


.39210315 


* 


102 


.90549867 


* 


1 


.12157650 


* 




.01089909 


* 


1 


. 17765968 


* 
* 


23 


.35792938 


* 


2. 


,12100000 


* 


7. 


,00000000 


* 







* 


MBrA : 


report on the main Weight 






* 


* 


90875 


no of raw events 


Ntot 


bl 


* 


* 


1000 


accepted events 


Nacc 


b2 


* 


* 





wt<0 events 


Nneg 


b3 


* 


* 


1 


wt>WTmax events 


Nove 


b4 


* 


* 


1.17765968 


WTsup, largest WT 


WTsup 


b5 


* 


* 


.00000195 


<Wt-WtMax> Overfl. 


AvOve 


b6 


* 


* 


.00000000 


<Wt> for Wt<0 


AvUnd 


b7 


* 


* 


.00018028 


AvOve/<Wt> , WT>WtMax 


ROverf 


b8 


* 


* 


.00000000 


AvUnd/<Wt>,Wt<0 


RUnder 


b9 


* 



*************************************************************************** 



MBrA: 



Detailed statistics for all branches 



KF 


AveWt 


ERela 


WtSup 


Wt<0 


Wt>Wmax 


Ntot 


Nacc 


Nneg 


Nove 


Nzer 


1 


.017100 


.024298 


.7062 


.000000 


.000000 


10024 


181 








96 


2 


.005937 


.026491 


.6685 


.000000 


.000000 


29851 


184 








343 


3 


.024082 


.024922 


.7078 


.000000 


.000000 


6738 


157 








57 


4 


.034175 


.026378 


1.178 


.000000 


.001035 


5022 


149 





1 


53 
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5 .046847 .022320 .7777 .000000 .000000 3332 168 38 

13 .002429 .036648 .8380 .000000 .000000 32573 91 474 

15 .020322 .038018 .6182 .000000 .000000 3335 70 40 

All: .010844 .010899 1.178 .000000 .000180 90875 1000 1 1101 



skipping some lines 

-630 Tau Pair: wtl, Spin Imprint weight 

nent sum bmin bmax 

281 .00000E+00 .00000E+00 .23000E+02 

undf ovef sumw avex 

.00000E+00 .00000E+00 .28107E+03 .13707E+01 





.0000 


.900000D+01 


oxxxxxxxxxxxxxxxxxxxxxxxxx 


I 




.1000 


. 120000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 




.2000 


. 190000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 




.3000 


.110000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 




.4000 


.900000D+01 


oxxxxxxxxxxxxxxxxxxxxxxxxx 


I 




.5000 


. 230000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 




.6000 


. 150000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 




.7000 


. 140000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 




.8000 


. 160000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 




.9000 


. 230000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


1 


.0000 


.210000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 


1 


. 1000 


. 150000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 


1 


.2000 


. 140000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 


1 


.3000 


. 170000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 


1 


.4000 


.110000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 


1 


.5000 


. 100000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 


1 


.6000 


.600000D+01 


oxxxxxxxxxxxxxxxxx 


I 


1 


.7000 


. 100000D+02 


oxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


I 


1 


.8000 


.300000D+01 


oxxxxxxxx 


I 


1 


.9000 


.200000D+01 


oxxxxx 


I 


2 


.0000 


.500000D+01 


oxxxxxxxxxxxxxx 


I 


2 


.1000 


.700000D+01 


oxxxxxxxxxxxxxxxxxxxx 


I 


2 


.2000 


. 100000D+01 


oxx 


I 


2 


.3000 


. 100000D+01 


oxx 


I 


2 


.4000 


. 100000D+01 


oxx 


I 


2 


.5000 


.000000D+00 





I 


2 


.6000 


. 100000D+01 


oxx 


I 


2 


.7000 


.300000D+01 


oxxxxxxxx 


I 


2 


.8000 


.000000D+00 





I 


2 


.9000 


. 100000D+01 


oxx 


I 


3 


.0000 


.000000D+00 





I 


3 


.1000 


. 100000D+01 


oxx 


I 


3 


.2000 


.000000D+00 





I 


3 


.3000 


.000000D+00 





I 


3 


.4000 


.000000D+00 





I 


3 


.5000 


.000000D+00 





I 


3 


.6000 


.000000D+00 





I 


3 


.7000 


.000000D+00 





I 


3 


.8000 


.000000D+00 





I 


3 


.9000 


.000000D+00 





I 



***************************************************************** 



* Tau Pair Finalize * 

* 1.00026167 +- .03630240 Spin Imprint <wtl> wtlave al * 

* 3.16386575 Maximum value wtl wtlmax a2 * 
*************************************************************************** 
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